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We present a Kadanoff-Baym formalism to study time-dependent phenomena for systems of interacting elec¬ 
trons and phonons in the framework of many-body perturbation theory. The formalism takes correctly into 
account effects of the initial preparation of an equilibrium state, and allows for an explicit time-dependence 
of both the electronic and phononic degrees of freedom. The method is applied to investigate the charge 
neutral and non-neutral excitation spectra of a homogeneous, two-site, two-electron Holstein model. This is 
an extension of a previous study of the ground state properties in the Hartree (H), partially self-consistent 
Born (Gd) and fully self-consistent Born (GD) approximations published in Ref. 1. We show that choosing a 
homogeneous ground state solution leads to unstable dynamics for a sufficiently strong interaction, and that 
allowing a symmetry-broken state prevents this. The instability is caused by the bifurcation of the ground 
state and understood physically to be connected with the bipolaronic crossover of the exact system. This 
mean-field instability persists in the partially self-consistent Born approximation but is not found for the fully 
self-consistent Born approximation. By understanding the stability properties, we are able to study the linear 
response regime by calculating the density-density response function by time-propagation. This functions 
amounts to a solution of the Bethe-Salpeter equation with a sophisticated kernel. The results indicate that 
none of the approximations is able to describe the respone function during or beyond the bipolaronic crossover 
for the parameters investigated. In overall, we provide an extensive discussion on when the approximations 
are valid, and how they fail to describe the studied exact properties of the chosen model system. 
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I. INTRODUCTION 


Many-body perturbation theory is one of the most 
common methodologies used to study quantum trans¬ 
port problems in which interactions among charge car¬ 
riers or between them and other constituents play a sig¬ 
nificant role. The method is based on diagrammatic per¬ 
turbation theory for the non-equilibrium Green’s func¬ 
tions together with a set of standard approximations to 
describe the many-body effects 2 . Although these approx¬ 
imations have been widely used, and thus their proper¬ 
ties explored, in the case of steady-state transport ', 
much less is known on their performance in the explic¬ 
itly time-dependent case ~ . This is particularly true 
for systems with moderate to strong electron-phonon in¬ 
teractions in which interesting phenomena like bistabil¬ 
ity and hysteresis have been observed . As such phe¬ 
nomena are typically driven by many-body interactions, 
it is natural to ask whether or not the approximate 
method can describe the relevant physics even qualita¬ 
tively. This is the case, for example, with the aforemen¬ 
tioned bistability whose existence has been subject to 
doubt on the quality of the method itself “ . The re¬ 


cent efforts to realize more sophisticated, but computa¬ 
tionally more demanding, approximations have enabled 
addressing these questions also in the framework of time- 
dependent many-body perturbation theory “ . It is im¬ 
portant to study these approximations on a wide scope 
to understand when they are predictive and the results 
can be trusted. Time-dependent many-body perturba¬ 
tion theory has also been recently applied to study vi¬ 
brational effects in ab-initio charge carrier dynamics in 
semiconductors e.g. relaxation processes after a laser ex¬ 
citation “ . This has become possible as further sim¬ 
plifications, in particular the generalized Kadanoff-Baym 
Ansatz 2, (GKBA), have been developed to keep the ap¬ 
proach computationally tractable along with the grow¬ 
ing system sizes. One could also in this manner study 
time-dependent phenomena in realistic molecular sys¬ 
tems continuing along the lines of the early studies of 
vibrational effects in photoelectron spectra of molecules 
In this context, in order to understand the reasons behind 
the successes or failures of the methods, it is neccessary 
to understand the many-body approximations underly¬ 
ing the additional simplifications such as GKBA. There 
are also topical fields in optoelectronics, such as cavity 
quantum electrodynamics, and optomechanics in which 
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one deals with formally similar systems as in the quan¬ 
tum transport case. Time-dependent many-body per¬ 
turbation theory has been used in these fields e.g. to de¬ 
rive time-dependent density functionals with preliminary 
results giving an indication of their quality-’' 1 ’ . There 
is however even less known about the properties of the 
approximations than in the more established quantum 
transport setup. 

In this work, we present an extension of a previ¬ 
ously introduced numerical method • to study time- 
dependent, inhomogeneous systems of interacting elec¬ 
trons and phonons. This is also an extension of the 
equilibrium formalism which we introduced in our ear¬ 
lier work in Ref. . Our approach is a variant of 
time-dependent many-body perturbation theory based 
on the Kadanoff-Baym equations (KBE) ’ . Here we in¬ 
troduce the relevant equations, time-dependent many- 
body approximations, and discuss some of their char¬ 
acteristic features e.g. the mean-field, Hartree (H) ap¬ 
proximation is shown to lead to the semi-classical Ehren- 
fest equations. The time-dependent partially (Gd) and 
fully- (GD) self-consistent Born approximations are in¬ 
troduced to study correlation effects beyond the mean- 
field level. These many-body approximations are in par¬ 
ticular suited to study time-dependent quantum trans¬ 
port with electron-phonon interactions as they are par¬ 
ticle number conserving in the sense of Baym ’In 
the present work, the method is instead applied to a fi¬ 
nite system since this allows us to assess its performance 
by comparing the approximate results to an exact solu¬ 
tion. Although the method can handle complex time- 
dependent perturbations, we restrict ourselves here to 
linear response functions obtained by tinre-propagating 
the Kadanoff-Baym equations. The density response 
function Sn/Sv obtained in this manner is equivalent to 
a solution of the Bethe-Salpeter equation (BSE) with an 
integral kernel which is a functional derivative 5Y*/5G 
of the self-energy E with respect to the electron prop¬ 
agator G '' ~ .The kernel therefore consists of dressed 
propagators and is fully frequency-dependent in the Born 
approximations. This level of approximation has, to our 
best knowledge, not yet been reached in the standard 
frequency-domain approach even for the much studied 
purely electronic systems “ . Moreover, as the phonon 
propagator is determined by the electronic response func¬ 
tion, by comparing to the equilibrium phonon propaga¬ 
tor we are able to comment on whether or not this ad¬ 
ditional level of sophistication amounts to improved re¬ 
sults. Lastly, we would like to note that although this pa¬ 
per is geared towards electrons and phonons, the method 
is in fact applicable to a variety of systems of inter¬ 
acting fermions and bosons e.g. electron-photon models 
(Rabi ,45 ) and electron-plasmon models (Lundqvist 16 ). 

As the application, we study a homogeneous, two-site, 
two-electron Holstein model which is a standard model 
describing interacting electrons and phonons 1 . This 
continues along the lines of our prior work in Ref. 1 
in which we focused on ground-state properties and 


studying the localizing effect of the electron-phonon in¬ 
teraction by comparing the many-body approximations 
against exact benchmark results. As a result, we found 
that the self-energy approximations gave rise to spon¬ 
taneous symmetry-breaking characterizable by an asym¬ 
metric electron density and nuclear displacement. The 
symmetry-broken solutions as well as solutions obtained 
by enforcing symmetry were analyzed with the help of 
total energies, energy components, and natural occu¬ 
pation numbers. It was concluded that the symmetry¬ 
breaking can be seen physically to mimic the bipolaronic 
crossover of the underlying system in which two nearly 
free electrons form a bound pair with an accompanying 
nuclear displacement. Moreover, out of the symmetric so¬ 
lutions, only the fully self-consistent Born approximation 
showed evidence of partially describing this crossover. 
Here we instead investigate the equilibrium electron and 
phonon propagators, and linear response functions of the 
same system using time-dependent many-body pertur¬ 
bation theory. The equilibrium propagators are studied 
in frequency-domain which gives a more detailed view to 
the properties of the approximations, and allows us to re¬ 
evaluate the physical picture obtained from the various 
energies. The linear response calculations on the other 
hand allow us to understand better the nature of the 
symmetric and asymmetric solutions found in our earlier 
work. In particular, we show that they are equivalent 
to the equilibrium solutions of the semi-classical equa¬ 
tions of the Dicke model 1 in which the appearance of the 
asymmetric solution represents the super-radiant phase 
transition in the thermodynamic limit “ . This transi¬ 
tion moreover appears as a bifurcation leading to insta¬ 
bility of the symmetric solution which in a finite system 
is not in agreement with the exact solution. One of the 
open questions addressed in this work concerns the stabil¬ 
ity of the symmetric and asymmetric solutions when go¬ 
ing beyond the mean-field approximation. In particular, 
we answer to the question whether or not the symmetric 
solution retains its stability in the Born approximations. 
Once stable solutions have been identified, we turn our 
attention to the linear response functions which are used 
to assess how the approximations describe the system 
reacting to a weak perturbation. There is a lot of sys¬ 
tematic work on static i.e. zero-frequency susceptibilities 
of either finite clusters 2 , or extended finite 53 and infi¬ 
nite 1 ’ dimensional systems with a focus on e.g charge- 
density wave phase transition temperatures. Here we 
thus extend these studies beyond the static case by con¬ 
sidering fully frequency-dependent response functions of 
a finite system. 

The paper is organized as follows. In Sec. II, we intro¬ 
duce our method: time-dependent many-body perturba¬ 
tion for electrons interacting with phonons. The method 
is applied in this work to the model system introduced 
in Sec. III. The results containing both the equilibrium 
electron addition and removal, as well as neutral exci¬ 
tation spectra are presented, analyzed and discussed in 
Sec. IV. The conclusions and an outlook are given in 
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Sec. V, and some more technical details are presented in 
App. A and B. 

II. THEORY 
A. Hamiltonian 

In the present work, we introduce the central concepts 
of time-dependent many-body perturbation theory for 
systems of electrons interacting with phonons. Although 
we do not discuss here electron-electron interactions, they 
could be included without additional conceptual difficul¬ 
ties. The time-dependent Hamiltonian operator is then 
given by 

&{z) + £ (fp( z )°l + fp( z ) & p) 

P P 

+ £M*)*fo 

ij 

+ £ £ ( m ifc( 2 ) a p+ mP kj( z ) & p) d ] d k, 

ij p 

where the properties of the system are encoded in the 
phonon frequencies ui p , generalized forces / p , elements 
of the electron one-body Hamiltonian Ay , and electron- 
phonon interaction elements mA. 

These quantities all depend on a time-argument ly¬ 
ing on the extended Keldysh contour 2 shown in Fig. 1. 
The Hamiltonian operator is however the same on the 
forward (—) and backward (+) branches, and indepen¬ 
dent of the contour time on the vertical equilibrium (M) 
track as in our prior work in Ref. 1. An explicit time- 
dependence allows one to realize a variety of physical 
scenarios from electrons and nuclei driven by electro¬ 
magnetic fields to more abstract simulations based on 
an interaction quench. In the present work, we however 
focus on another type of time-dependence arising from 
the choice of the initial state. 

The electrons and phonons are described in second 
quantization with annihilation Cj, a p and creation cj, aj) 
operators which obey canonical anti-commutation and - 
commutation relations, respectively. In order to facilitate 
a compact presentation of the many-body perturbation 
theory, we further introduce the self-adjoint phonon op¬ 
erators 

01 ,p = (hjj T ^p) /, 02,p — (0» ^p) /, 

to which we associate a collective index P = € 

{1,2 },p} so that we can write their commutation rela¬ 
tion compactly as 

[0P,0Q] = OtpQ , 

where 0/.\ p \q &2p,2q 0 and Oi\ p2 q &2q,lp l&pq' 

These operators can be physically understood as com¬ 
ponents of the displacement (0i p ) and momentum (02 P ) 



FIG. 1. The extended Keldysh contour which consists of 
the vertical, imaginary-time track responsible for the initial 
equilibrium preparation, and of the horizontal forward (—) 
and backward (+) real-time tracks related to the real-time 
time-evolution, (color online) 


operators. They allow us to rewrite the Hamiltonian op¬ 
erator as 

£(*)=£ ^ pq (~) 0 p 0 q + £ Fp{z)(j)p 

PQ P 

+ £ h iji z )c\cj 
ij 

+ ££ M b( z )^ pS i^» w 

ij P 

where the phonon frequencies, generalized forces, and 
electron-phonon interaction are incorporated into 

-^ p ,2(/p(V) - fp(z))/V 2, 

^P<ip,q c i q ( z ) = w p(' 2: )(^pi?^p‘:, + a P‘J P ,9‘:,)/2 j 

M jk( z ) = KA mP jk( z )+ m Tj( z ))/^ 

-iS, pt2 (m^ k (z) -m p k *(z))/V 2 , 

which are to be understood in this work to represent el¬ 
ements of a vector, matrix, and a vector of matrices, re¬ 
spectively. The one-body electron Hamiltonian elements 
are also to be understood as elements of a matrix. In 
the following an overhead arrow denotes a vector (F), 
boldfaced symbols matrices ($ 7 , h ), and a combination of 
these two a vector of matrices (Af), while tr denotes a 
matrix trace. 


B. Many-Body Perturbation Theory 

The central quantities of many-body perturbation the¬ 
ory of interacting electrons and phonons are the phonon 
field expectation value, and the phonon and electron 
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propagators defined as 


(j>p{z) = — Tr 
Dpq{z;z') = ^Tr 


Gij(z-,z') = — Tr 


r | e -z/ c dSH(2)^ p ( z )j 


T { 

T { 


g ~ i Ic dz H ( Z 


g ~ i Ic dz H ( Z '> 


'A<t>p(z)A(j>Q(z ')}! , 


c^cjO')} 


( 2 ) 

(3) 


where A<)>p = (f>p — (f>p is a fluctuation operator, Z = 
Tr[e~ l J dz H ( z )j the partition function, Tr the trace over 
a complete set of quantum states, and T is the time¬ 
ordering operator on a Keldysh time-contour G of Fig. 1 
acting on operators given in the Schrodinger picture but 
having time-arguments z, z' for book-keeping reasons . 
These objects have a closed form perturbation expansion 
obtained using Wick’s theorem and re-summing all terms 
into two electron and phonon propagator line irreducible 
contributions. This leads to the equations 


$(z) = f dz d(z 7 z) (F(z) — itv{M(z)G{z\ z + )) , 

Jc 


(4a) 


D(z] z') = d(z; z') + I dzdz'd(z; z)Tl(z; z')D(z'; z r ), 

Jc 


(4b) 


G(z; z') = g(z 7 z!) + f dzdz' g(z 7 z)'S(z; z')G(z'; z!). 
Jc 


(4c) 


where g and d denote the non-interacting electron and 
phonon propagators defined by Eqs. (3) and (2) in the 
absence of the electron-phonon interaction. The inte¬ 
gral kernels E = £[G, D] and II = II[G, D\ are non-local 
one-body potentials known as electron and phonon self¬ 
energies. These self-energies contain information on in¬ 
teractions of the system, as well as the external driving 
induced by the generalized force F. The non-interacting 
electron and phonon propagators are given respectively 
by 


g{z\z') = -iU(z,t 0 )(o(z,z') - f + (/3h M )^)U(t 0 ,z'), 
d(z; z') = -ia.V{z,t 0 )[e{z,z') + f_(f3n M a)^V(t 0 ,z '), 

where 6 = 6 1 with 6 being the Heaviside function and 
1 the identity matrix, /3 is the inverse temperature, f± 
denote the Fermi-Dirac (+) and Bose-Einstein (—) dis- 

- M 

tribution functions, fl = fI(to — w) independent of r 
is the Matsubara component of 

J7(z) = $ 7 ( 2 ) + ft 1 (z). 


Finally, we introduced the time-evolution matrices as so¬ 
lutions to 


id z U(z,z') = h(z)U(z, z '), 

-idgi U(z, z') = U(z, z')h(z '), 
id z V(z, z') = Cl(z)aV(z, z '), 
-id z >V(z,z’) = V(z,z')£l(z)a, 

with the initial conditions U(to 7 to) = V(to 7 to) = 1. 

In our earlier work in Ref. 1, we introduced our imple¬ 
mentation of the equilibrium Matsubara formalism ob¬ 
tained by choosing time-arguments z = to — it, z' = 
to — it ' on the imaginary track. Here we focus on an 
extension of this formalism to time-dependent cases in 
which it is more natural to differentiate Eqs. (4) with re¬ 
spect to the first contour time in order to arrive at the 
equations of motion 

(■ iocd z - Cl(z))$(z) 

= F(z)-itr(M(z)G(z-,z+)), (5) 

( iad z — £l(z))D{z- 7 z') 

= S(z,z')+ f dz II(z; z)D(z; z') , (6) 

Jc 

( id z - h(z))G(z\ z') 

= 8(z, z') + f dz 53 ( 2 :; z)G(z; z'), (7) 

Jc 

where i = 1 1. These equations together with their con¬ 
jugate equations obtained by differentiating with respect 
to the second time-argument of the propagators form a 
closed set of the equations which can be solved once an 
approximation for the many-body part of the self-energy 
has been fixed. 


C. Self-Energies 

The self-energy E, as noted above, contains both a 
contribution arising from the generalized force Fp(z), as 
well as a part induced by the electron-phonon interac¬ 
tions. The phonon propagator, being defined in terms of 
fluctuation operators, is not directly influenced by this 
force, instead it appears in the electron self-energy and 
can be handled by writing the self-energy as 

£«(*;*') = 6(z,z')v ni ij(z) + Emb ,ij(z;z r ) 
where we introduced the potential 



which represents the classical potential induced by nu¬ 
clei experiencing a generalized force Fq . The many-body 
self-energy, denoted by MB, is then subject to approxi¬ 
mation. The approximations used here, and introduced 
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earlier in Ref. 1, are summarized diagrammatically in 
Fig. 2. The approximate electron self-energies consists 
of the Hartree (H) and Fock (F) diagrams. The Hartree 
diagram can be written as 

E H [G\ij{z;z') = S(z, z')vn[G]ij{z), 

where the time-local Hartree potential is given by 

vh[GUz) = -iY^Mij(z) 

kl 

PQ 

x [ dz d P Q (z, z)M® (z)Gik (z;z + ). (8) 

Jc 

or alternatively by 

vn[G\ij{z) = E W<Ma) - v njij (z) , (9) 

p 

which follows from the equation of motion for the non¬ 
interacting phonon propagator. Electron self-energy 
terms beyond Hartree contribute to the exchange- 
correlation, many-body self-energy 

E X c ,ij(z’i Z ) = ^MB Z ) ^H,ij (zj Z ) j 

whose lowest order diagram is the Fock diagram 

E F [G, £>]#(*:; *0 

= * E Mg{z)M%{z')Dp Q (z-,z')G kl {z-,z'), 

kl,PQ 

which is a time-nonlocal memory term describing single¬ 
phonon absorption/emission processes. The only phonon 
self-energy diagram used in this work is the bubble dia¬ 
gram 

n B [G]pQ(£; z') 

= -*E ( z ) M ki( z ') G ii{ z '-, z)G jk (z ; z'), 

ijjkl 

which describes simple phonon induced electron-hole 
excitation processes. 

The many-body self-energies and their abbreviations 
used throughout the text are summarized in the list be¬ 
low. 

H: The Hartree approximation consists of approximating 
the electron self-energy with the Hartree diagram 

E H [G\ij(z m , z') = S(z, z')vn[G\ij(z ), 

and neglecting the phonon self-energy. This is a 
mean-field approximation in which electrons feel 
only the classical potential due to nuclei. The re¬ 
sulting Hartree equations 

i^-G{z-,z + ) = [h{z) + v n (z) 


+ vn(z),G(z-,z + )\ , 

(10a) 

iad z $(z) = Cl(z)(j)(z) + F(z) 


— itr(MG(z; 2 i + )) , 

(10b) 



FIG. 2. The Hartree (H), partially self-consistent (Gd), and 
fully self-consistent (GD) Born self-energies summarize the 
many-body approximations used in this work. A two-fold line 
with an arrow indicates a dressed electron propagator, while 
single and two-fold wiggly lines represent bare and dressed 
phonon propagators, respectively. An open circle and a closed 
circle represent a connection for a phonon and an electron 
propagator, respectively. 


can be shown to be equivalent to the semi-classical 
Ehrenfest equations, see App. A. 

Gd: The partially self-consistent Born approximation 
amounts to approximating the electron many-body 
self-energy with 

E Gd [G\ij(z]z') = Tj-a[G\ij(z\ z') 

+ Ep[G, d\ij(z; z '), 

where d is the bare phonon propagator obtained 
by putting the phonon self-energy to zero. This 
amounts to saying that that the nuclei are unaf¬ 
fected by the electronic particle-hole excitations. 

GD: The fully self-consistent Born approximation is de¬ 
fined by writing the electron many-body self-energy 
as 

£ gd [G, D] iJ (z-,z') = E h [ Gjij&z') 

+ Ep[G, D\ij(z\ z 1 ) 

while the phonon self-energy is given by 

n G D [G]pq(z; z ') = n B [G]pQ(z; z l ). 

Note that although we dress the phonon propagator 
in the Fock diagram, one should not use a dressed 
propagator in the Hartree diagram as it leads to 
double-counting of terms in the perturbation ex- 

• 9 

pansion. 
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- $Gd = 




-$gd = 




FIG. 3. The ‘F-functionals for the Hartree (H), partially self- 
consistent (Gd), and fully self-consistent (GD) Born approx¬ 
imations. A two-fold line with an arrow indicates a dressed 
electron propagator, while single and two-fold wiggly lines 
represent bare and dressed phonon propagators, respectively. 
Note that the minus sign on the left hand side arises due to 
the loop rule . 


These approximations are all $-derivable, that is the cor¬ 
responding self-energies can be obtained as the functional 
derivatives 


X lJ {G,D](z-,z') = 


6$[G,D] 


SGji(z';z ) ’ 


U PQ [G,D](z]zf) = -2 


5$[G,D} 


SD qp (z';z) 


of an approximate ^-functional which are shown in 
Fig. 3. Note that the subscript S refers to a symmetrized 
functional derivative, see . The <F-derivability of the elec¬ 
tron self-energy together with self-consistency in the elec¬ 
tron propagator guarantee gauge invariance and conse¬ 
quently fulfillment of the electron density conservation 

law 35 - 36 . 


D. Kadanoff-Baym Equations 


The equations of motion of Eqs. (5) are customary 
solved by projecting the propagators to different parts 
of the Keldysh contour by choosing the time-arguments 
appropriately . This procedure leads to the greater (>), 
lesser (<), left ([), right (]), and Matsubara (M) com¬ 


ponents 


a^{t;t’) = a{t±-t ' T ), 
a) (t; t) = a(t; to — it) , 
a^(r; t) = a(to — it\ t ), 
a M (r; t') = a(t 0 - it ; t 0 - it') , 

where the subscript y denotes a time evaluated on the 
forward/backward branches of the contour, and a is a 
function in the space of Keldysh functions 3 . The Keldysh 
components of the phonon and electron propagator obey 
the symmetries 

Gf j (t,t') = -[Gf i (t';t)}*, 

G\j{T-,t) = [G]i(t-, 0 - t)]* , 

D%W) = -[tfJpCGC = -[Df Q (t,t')]\ 
d pq(t; t) = [D ] Qp {t; (3 - t)]* = [Dp Q (/3 - r; t)] * , 

where the additional symmetries of the phonon propaga¬ 
tor are due to the symmetry Dpq(z;z') = Dqp(z! , z). 
The equations of motion obtained by taking all possible 
components of the contour-time equations of motion form 
a set of non-linear integro-differential equations of motion 
known as the Kadanoff-Baym equations . The symme¬ 
tries of the propagator however imply that we only need 
equations of motion for the greater, lesser, and right com¬ 
ponents where the first two are required for times t >t r . 
The relevant equations of motion are then 

idtG% (t; t') = h cff (t)G * (t; 0 + P? [£ xc , G] (i; t') , 
id t G\ (t; r) = h oS (t)G 1 (i; r) + f [£ xc , G](f; t') , 
id t D * (t; t') = ol ( p(t)D * (i; t') + I* [n, D] (t; t')) , 

id t D\f,T) = a(Cl(t)D\t-,T) + l\n,D]{f,T)) , 
for off time-diagonal and 

*^G*(i;i) = h eS (t)G*(t;t) + I^[£ xc ,G](f;f) + h.c., 

ij t D^{t-,t) = a(Cl(t)D^(t;t) + I^[H,D](t;t)^j + h.c., 

where h.c. denotes the Hermitian conjugate, for on time- 
diagonal time-propagation. Here we introduced the ef¬ 
fective one-body electron Hamiltonian 

h c s(t) = h(t) + v n (t) + v H (t), 

as well as the collision integrals 

I^[a,b](t;t') = [o' *tf](t;t') 

+ [a R • b^] (t; t') + [a^ • b A ] (t; t'), 

I 1 ' [a, 6](t; r) = [o' * b M ] (t; r) + [a R • b 1 ] (f; r). 
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with the bullets and stars denoting convolution integrals 
of the form 

/*oo 

\a • b\(t\t') = / dt a(t,t)b(i,t '), 

Jt 0 

rP 

[a * b\(t;t') =—i / dr a(t,T)b(T,t'), 

Jo 

where a and b are possibly matrix valued functions on 
the Keldysh contour. The Hartree potential appearing in 
the effective Hamiltonian can be evaluated using Eq. (9) 
instead of Eq. (8) by taking advantage of the equation of 
motion 

idt$(t) = a(Cl(t)$(t) 

+ F(t) - , 

for the real-time <j>p(t ) = <fip(t±) phonon field expec¬ 
tation value. The Kadanoff-Baym equations, including 
the equation above, then form a closed set of equations 
which, when supplemented with the initial conditions 

Gf j (t 0 -,t 0 ) = G%(0 ± ), 

G ] ij{to-,T ) =G.f/(-r), 

Dp Q (to;t 0 ) = DpQ ^), 
D], Q (to-,T)=D? Q (-T), 

<t>p(to) = J>p , 

given by the equilibrium Matsubara components intro¬ 
duced in Ref. 1, can be solved on a computer by time 
propagation . 

III. MODEL 

Our model system is a two-site Holstein model 
which can be viewed as a minimal representation of a 
system in which electrons move between two molecules, 
so that they are coupled to the local vibrational modes of 
these molecules. In the case of two identical molecules, 
we find that only the relative displacement couples to the 
electron density difference between the molecules, and 
thus the Hamiltonian operator for the isolated system 
reduces to 

H m = wo aJa 

~ ^kin ^ ' (ci a C2a + C 2cr Cicr) 
a 

- -^=(a f + a) ^2 (™io- - « 2 a) , 

where a and a' annihilate and create a phonon to the 
relative displacement mode, Ci a is the electronic opera¬ 
tor that annihilates an electron of spin a at site i, and 
hie = c\ a c i(7 is the electron density operator at site i. 


The parameters wo, ikin and g characterize the bare vi¬ 
brational frequency, inter-site hopping and local electron- 
phonon interaction strength, respectively. This Hamilto¬ 
nian corresponds to the relative Hamiltonian of Ref. 
which is chosen here over the full Hamiltonian due to 
both its simplicity and computational reasons. The sys¬ 
tem can be probed with external time-dependent fields 
which are described with the Hamiltonian operator 

H(t) = H m + f{t){a) + a) + ^2,Vi{t)h ia , 

icy 

where / and v. t describe amplitudes of the external fields 
acting on the nuclei and electrons, respectively. The 
displacement and momentum operators, defined in this 
model as u = (fit + a)/y/2 and p = «(fit — a)/y/2, allow 
us to rewrite the Hamiltonian operator as 

H(t) =H M + V2f(t)u + YMt)n t g, (Ha) 

icr 

^kin ^ ^ (^lcr^2cr “I” 

(T 

- gu^2 (hia - h 2 a), ( 11 b) 

a 

which is equivalent to the Hamiltonian of Eq. (1) with 
the matrix elements 

K P {z) = 6(t 0+ ,z)S, Z) iV2f{z ), 

?<,(*) = ( Ks q + a ^ q )/ 2 > 

“ 5 a) > 

where we dropped the phonon mode index due to hav¬ 
ing only the relative mode. Moreover, here 9 denotes 
a contour-time Heaviside function and the origin of 
the backward branch. The time-independent properties 
of this model depend on the two dimensionless parame¬ 
ters 

w 0 

7= —, 

^kin 

t'kin^O 

denoting the adiabatic ratio and effective interaction. 
The adiabatic ratio 7 describes the relative energy scale 
of electrons and nuclei, while the effective interaction A is 
a measure of the coupling between the motions of these 
two constituents. 


IV. RESULTS 

In the following, we present our results for the equi¬ 
librium propagators and linear response functions. The 
results are for a system initially at zero temperature in 
the pure two-electron N = 2 spin singlet S 2 = S z = 0 



ground state. This is mimicked in many-body perturba¬ 
tion theory with the inverse temperature = 10 3 . 

Moreover by choosing G ia ^ a \z\ z') = 5 arT 'Gij(z ; z') such 
that N = —2 tJ2iG^(t;t) = 2 for all times, we can en¬ 
sure that S z = 0. The results cover the physical pa¬ 
rameters 7 = 1/2,1/4 and A G [0,2] corresponding to 
the weak- and intermediate-to-strong interactions. The 
approximate results (H, Gd, GD) are obtained by first 
solving the imaginary-time Matsubara propagators with 
an imaginary-time grid, solution method and related pa¬ 
rameters identical to the ones used in our previous work 
in Ref. 1. This leads to multiple solutions character¬ 
ized by either symmetric or symmetry-broken electron 
densities and nuclear displacements as shown in Ref. 1. 
The former kind are the only solutions for a sufficiently 
weak interaction and are known here as the symmet¬ 
ric solutions, while the latter kind arise for sufficiently 
strong interactions and are referred to as asymmetric so¬ 
lutions. Here we mention that our approximations do not 
respect the exact transformation relating the relative and 
full Hamiltonians of Ref. . This is seen as quantitative 
differences between some equilibirum observables, which 
are invariant under this transformation in the exact case, 
presented here and in Ref. 1. In the present work, the 
real-time electron and phonon propagators are then cal¬ 
culated by time-propagating the Kadanoff-Baym equa¬ 
tions, according to an adapted version of the algorithm 12 , 
using the abovementioned Matsubara propagators either 
directly or indirectly (see the linear response section) 
as initial values. The time-grid is uniform with a grid 
spacing or time-step St such that fkin^T G [0.025, 0.075] 
extending from zero to the final time T chosen so that 
tk; n T = 200. The time-domain propagators are finally 
Fourier transformed to arrive at their frequency-domain 
representations. The Fourier transforms are calculated 
with a high-order quadrature formula and unless other¬ 
wise stated by using the Hanning window function . 


A. Equilibrium Propagators 


The out-of-equilibrium behavior of a system can be 
better understood if we first understand the equilibrium 
properties of this system. These properties are deter¬ 
mined by the equilibrium electron and phonon propaga¬ 
tors which we have studied in Ref. from the perspec¬ 
tive of time-local (e.g. density matrix) and integrated-out 
(e.g. total energy) observables. Here, we further shed 
light on the quality of our approximations by investigat¬ 
ing the frequency structure of these propagators. In this 
section, the propagators depend only on the relative time 
and our convention for evaluating Fourier transforms is 
that the first time argument is integrated over and the 
second kept fixed at the initial time. 


1. Electron Propagator 


The electron propagator is directly related to the pho- 
toemission i.e. electron removal and inverse photoemis¬ 
sion i.e. electron addition spectra. This can be qualita¬ 
tively seen from its zero-temperature Lehmann represen¬ 
tation 




where n(/ ±:L = E^ ±l — Eq is the electron addi¬ 

tion/removal energy while = (\k)/ +1 | c] CT I’I'oO and 
fnia — _1 | Gct |'E'o r ) are the corresponding ampli¬ 

tudes. Here and E% denote the nth eigenstate 
and -energy of the N electron system. The Lehmann 
form is used below to interpret the results shown in 
Fig. 4 for exact diagonalization (ED) and many-body 
perturbation theory (H, Gd, GD). The greater and lesser 
components are related by the particle-hole symmetry 
Gfj(uj) = — ( — 1 ) l-J G^(—w) whose fulfillment is dis¬ 
cussed below, and therefore we only show results for the 
lesser component. Let us focus first on the main panels 
(contour plots) to illustrate the overall frequency struc¬ 
ture, and start by examining the exact results. The exact 
spectra develop as a function of the interaction from the 
singly peaked, non-interacting spectra described by the 
function 

= T«7r(Tl) i_i <5(a; =F i kin ), 

into spectra consisting of multiple peaks whose positions 
are up to an energy shift given by the energies of the 
one-electron system. The one-electron energies E^ =1 are 
nearly uniformly separated by the bare phonon frequency 
for a weak interaction A « 1. This manifests itself in 
the exact spectra as emergence of the so-called phonon 
sideband structure which gains intensity as the initial 
distribution loses intensity. In the case of a sufficiently 
strong interactions the lowest energies instead consist of 
nearly degenerate pairs separated by the bare phonon 
frequency ' 1 ’ . In this case the one-elecron system can be 
characterized as polaronic and is, as a first approximation 
in the limit A> 7 , described by 

|^ F ±CT ) = ^(cLl±cLl t )|0;fc) 


where X = exp(— igp/uio) is a shift operator, and | 0 ; k) 
is an empty electronic state and fcth eigenstate of a^a}' 2 . 
In the same limit, we find the two-electron ground state 

!^o F) = ^(4 t 4* 2 + 4 t 4* t2 ) |0;0) . 

which has a bipolaronic character. The removal energies 
and associated amplitudes 


&k,lcr ~ ^ kin A/4 + Ul 0 k , 


y, y, [fkda'y 

Z6{±} 


2 


e —A/4 7 / X \ 

2k\ \47y 
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then show that the spectra consist of peaks separated by 
the bare phonon frequency with intensities following a 
Poisson distribution . The exact results shown in Fig. 4 
indicate that the initial spectra become denser as inter¬ 
action is increased such that the two lowest excitations 
approach one another faster than the third which stays 
roughly a bare phonon frequency apart, especially for 
7 = 1/4. At the same time spectral weight is redis¬ 
tributed in particular to the third and higher-lying exci¬ 
tations. We interpret this as a precursor of the crossover 
to a Poissonian disribution which is a signature of a pola- 
ronic one-electron and bipolaronic two-electron system. 
This change is accompanied by an overall shift of the 
spectra to higher energies which appears smoothly as a 
function of the interaction, although more rapidly around 
A ~ 1 for the adiabatic ratio 7 = 1/4. The shift im¬ 
plies that one needs more energy to either add or remove 
electrons indicating that the two-electron ground state 
becomes more stable. This is in agreement with the in¬ 
crease in the bipolaron binding energy, see e.g. Ref 1, and 
is hence associated with the fact that the two-electron 
ground state becomes characterizable as bipolaronic. In 
addition to these changes there is a faint signal around 
aj/fkin ~ —3 for 7 = 1/2 and weak interactions, which 
is to be understood as the removal energy associated 
with the anti-bonding state of the one-particle system. 
This feature is washed out for the lower adiabatic ratio 
7 = 1/4, in contrast to a similar feature of the single¬ 
electron case" . 

The question is then how well the many-body approxi¬ 
mations reproduce the qualitative features of these spec¬ 
tra and thus the associated physics. The Hartree ap¬ 
proximation leads to spectra with peaks located at the 
eigenvalues of the equilibrium Hartree equations. In the 
case of the symmetric solution, the Hartree potential van¬ 
ishes, and this approximation just reproduces the non¬ 
interacting result gfj(co) for all values of the interaction 
thus failing to describe the exact spectra. The asymmet¬ 
ric case however displays a more complicated behavior 
with the propagators given by 

Gjj a+ ij( ljJ ) — *7rA S(uJ ^ tkinA) , 1 / j, 

G H a+ i*M = T*tt(i ± (-l)Vl - A~ 2 )<5(u; =F uA), 

where H a + denotes the asymmetric solution with a pos¬ 
itive density difference n\ a — ri 2 o- The asymmetric 
spectra emerge at A = 1, and contain a single peak 
which moves to higher energy linearly as a function of 
the interaction. The particle-hole symmetry is broken 
along with the reflection symmetry, however they are 
replaced by G^ a+ij (u) = G^ aji (u) and G£ o+ii (w) = 
— ( — 1 Y~-’G^ a _j i (—uj) where H a _ denotes the asymmet¬ 
ric solution with a negative relative density. These rela¬ 
tions represent the original symmetries under the inter¬ 
change of the two degenerate asymmetric solutions. Al¬ 
though the asymmetric solutions lead to spectra which 
shifts to higher energies as the exact spectra do. they do 


not show signs of the phonon sideband structure. The 
Born approximations correct this flaw and show a clear 
sideband structure. The partially self-consistent Born 
approximation, in the case of the symmetric solution, 
however shows that all removal energies behave roughly 
in a similar fashion, namely they increase monotonously 
and nearly linearly as a function of the interaction. The 
spectra do not show signs of a peak corresponding to a 
removal energy associated with the anti-bonding state of 
the one-particle system. This feature instead emerges 
qualitatively correctly in the fully self-consistent approx¬ 
imation. The fully self-consistent approximation also im¬ 
proves the position of the dominant removal energies for 
weak interactions by showing a stronger increase of the 
lowest removal energy and a simultaneously decrease in 
the sideband removal energies. Moreover, on the con¬ 
trary to the monotonous behavior of the partially self- 
consistent approximation, the fully self-consistent ap¬ 
proximation shows a signature of a stronger change in 
the structure of the spectrum for A = 1/4 approximately 
where the exact spectrum also changes. At this point, 
the fully self-consistent spectrum however becomes too 
dense, and does not shift correctly to higher energies. 
The asymmetric solutions, once they appear for a suffi¬ 
ciently strong interaction, are similar in these approxima¬ 
tions and differ from the asymmetric mean-field solution 
by the fact that there is a related sideband structure. 

The changes in the spectra from a non-interacting to 
a fully interacting case should emerge in a way which 
respects the two lowest order sum rules for the electron 
propagator 


OO 



— OO 


= N , ( 12 a) 

OO 

g (1) ^ / ^wtrG<( W ) 

—OO 

= E e + E ep , (12b) 

where the right-hand sides are equilibrium expectation 
values of the electron number N, and electron E e and 
electron-phonon interaction E ep energies, see Ref. 1. The 
top panels of Figs. 4 show that both constraints are ful¬ 
filled up to a numerical accuracy. The numerical devi¬ 
ations especially for 7 = 1/2 are due to choice of time 
discretization and frequency integration. Moreover, we 
note that all frequency moments have been calculated in 
the present work from spectra obtained using a rectangu¬ 
lar window function. The first moments, which are equal 
to the mean of the distribution, show that in addition 
to the asymmetric cases, only the exact and symmetric 
fully self-consistent spectra move appreciably to higher 
energies, in particular for 7 = 1/4. 

The left panels of Fig. 4 show the position and intensity 
of the lowest lying peak labeled with [E], as in Electronic, 
of the removal spectra. This peak is the most significant 
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FIG. 4. The exact (ED) and approximate (H, Gd, GD) electron propagator as a function of the interaction A and frequency 
uj. The top and bottom figures correspond to the adiabatic ratios 7 = 1/2 and 7 = 1/4, respectively. The main panels in the 
middle show |trG < (o))| /T for the symmetric (s) solutions and, starting from the critical interaction denoted with a vertical 
dashed line and its value Ac, for the asymmetric (a) solutions. The top panels show the zeroth g 1 ' 0 ' 1 and first g ^ moments 
(points) of the spectra and the corresponding expectation values (lines) to illustrate the fulfillment of Eqs. (12). The left panels 
show the intensities and positions of the lowest energy peak of 2 |GJj i (w)| /T and 2 |G^ 2 (^)| /T, labeled as [E] in the main 
panels, for the symmetric (solid line) and asymmetric (dashed line) solutions. The right panel show |trG < (oi)| /T on a linear 
scale for the symmetric solutions at A = 2. (color online) 


part of the spectra in the regime of weak to intermediate 
interactions where the many-body approximations are 
expected to be in qualitative, or even quantitative, agree¬ 
ment with the exact solution. Our results show that, out 
of the symmetric solutions, the Born approximations are 
indeed in a good agreement with exact results in the weak 
coupling regime. The partially self-consistent approxima¬ 
tion however deviates considerably already for interme¬ 
diate interactions A ~ 1, while the fully self-consistent 
approximation gives a reasonably good estimate of both 
the position and intensity up to borderline strong inter¬ 
actions A ~ 1.5. For stronger interactions, both approxi¬ 
mations fail to describe the shift of the position, as well 


as the decrease of the intensity correctly, although the 
fully self-consistent approximation gets the latter trend 
considerably better. The exact position and intensity of 
this peak change more abruptly in the case of 7 = 1/4, 
and imply that the sidebands become the most intense 
part of the exact spectra for the higher interactions con¬ 
sidered in this work. The many-body approximations do 
not show sufficient loss of intensity, and therefore fail to 
redistribute the spectral weight correctly to the higher 
energy part. The results for A = 2 shown in the right 
panels of Fig. 4 verify this statement and moreover show 
that the approximate spectra do not bear resemblance to 
the shape of the exact spectra. Lastly, the asymmetric 
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solutions capture the loss of the intensity qualitatively 
correctly for the site with the lower occupation but in 
doing so break the reflection symmetry which leads to 
an increase of the intensity of the site with a higher oc¬ 
cupation. This is natural since it becomes favorable to 
remove electrons from an already almost fully occupied 
site and vice versa. 

To summarize, we found that for the adiabatic ratios 
considered here the Hartree, and partially and fully self- 
consistent Born approximations are in a good agreement 
with exact results for very weak A -C 1, weak A < 1, 
and intermediate A ~ 1 interactions, respectively. More¬ 
over, the agreement between exact and approximate re¬ 
sults improves when the electronic and phononic energy 
scale become closer to one another for 7 = 1/2. These 
observations are similar to the conclusions of our ear¬ 
lier work in Ref. in which it was further observed that 
when approaching the anti-adiabatic limit the approxi¬ 
mate results start to again deviate from the exact results. 
In particular, the comparison of the total energies and 
natural occupation numbers conducted in our previous 
work supported the view that the fully self-consistent ap¬ 
proximation describes the bipolaron crossover partially. 
The present results show that as the interaction A is in¬ 
creased none of the approximate removal spectra i) move 
to higher energies as in ~ 3tki n A/4 nor ii) develop towards 
a uniformly uto-spaced distribution with a Poissonian-like 
envelope. The results are consistent with our earlier find¬ 
ings as the sum rules are satisfied and e.g. E e + E ep 
does show a clear significant increase in the fully self- 
consistent approximation. As discussed above, points 
i) and ii) signal a bipolaronic system, and their incor¬ 
rect description rather suggest the conclusion that none 
of the approximations describe the bipolaronic crossover 
even partially. The failure to describe ii) is related to the 
observation that the intensity of the lowest excitation 
energy does not decay fast enough as a function of the 
interaction in the approximate results. This is analogous 
to the insufficiently fast decaying quasi-particle spectral 
weight used as an indicator of absence of the bipolaronic 
metal-insulator transition in the fully self-consistent ap¬ 
proximation . Finally, our observation on the relation 
between the frequency-resolved and integrated-out quan¬ 
tities is similar to those obtained earlier e.g. for the GW 
approximation in the homogeneous electron gas in which 
self-consistent total energies were good but the plasmon 
description inadequate . 


2. Phonon Propagator 

The phonon propagator is an indicator of the proper¬ 
ties of the nuclear system, and relates to neutral excita¬ 
tions, as shown by its zero-tenrperature Lehmann repre¬ 


sentation 

^pq( w ) = Dq P (—u) 

= -*27T J2fn,pfn,Q*S(“-^) 
n 

where = E^ — Eq is a neutral excitation energy, and 
f^ P = (y%\A<j>p I'I’o) the corresponding amplitude. 
The frequency-domain phonon propagators obtained by 
means of exact diagonalization (ED) and many-body the¬ 
ory (H, Gd, GD) are shown in Fig. 5. Let us first discuss 
the contour plots which illustrate the overall frequency 
structure of the spectra. The exact results show that 
as the interaction is increased the initial, non-interacting 
spectra described by 

dpp(w) = —i2tt5(uj — wo), 

dp Q (u) = -2n(P - Q)6(w - w 0 ) , P , 

where P, Q £ {1,2} with 1 and 2 referring to the rel¬ 
ative displacement and momentum, develop into multi- 
peaked spectra consisting of a low and a high energy 
scale. The low energy scale consists for a sufficiently 
strong interaction of a single high intensity peak accom¬ 
panied by a weaker peak separated approximately by the 
bare phonon frequency. The high intensity peak which 
is labeled with [P] referring to Polaronic in the figures, 
develops continuously from the initial distribution and 
moves rapidly towards zero energy as a function of the 
interaction strength. This is true for both adiabatic ra¬ 
tios with the difference that [P] approaches zero more 
abruptly for 7 = 1/4. The high energy scale, on the other 
hand, consists of multiple low intensity peaks above the 
first electronic excitation energy of the non-interacting 
system. As the interaction is increased, these excita¬ 
tions move towards higher energies and, although ini¬ 
tially gain intensity, become suppressed for a sufficiently 
strong interaction. These features can be understood 
from the adiabatic potential energy surfaces defined and 
analyzed in Ref. 1 and shown here in Fig. 6 . This fig¬ 
ure shows that the initially quadratic lowest potential 
energy surface Eq(u) becomes more shallow as the inter¬ 
action is increased which is seen in Fig. 5 as a decreasing 
phonon frequency. The surface builds up a double-well 
structure for A > 1 which manifests itself in the exact 
results as a nearly degenerate ground and first excited 
state [P]. Moreover, as the barrier between the wells in¬ 
creases, the low energy spectra approach the harmonic 
spectra of the isolated wells which is seen in the exact 
results for A = 2 as a single peak located roughly at 
the bare phonon frequency. The high-energy spectra, on 
the other hand, agree with the first excited state surface 
Ei ( u ) remaining roughly quadratic while the surface sep¬ 
aration Pi (it) — E 0 (u) increases. As discussed in Ref. 1 , 
in the adiabatic case 7 < 1 the double-well structure is 
correlated with a splitting of the nuclear ground state 
probability distribution and the crossover to a bipola¬ 
ronic state. In this section, we thus identify its spectral 
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signature, that is the low energy part consisting of the 
two peaks, as an indicator of a bipolaronic state. 

Let us then focus on the approximate results. The 
Hartree and partially self-consistent Born approxima¬ 
tions approximate the phonon propagator with the non¬ 
interacting propagator which does not describe the true 
behavior of the interacting system discussed above. The 
question is then how the fully self-consistent approxima¬ 
tion, in which the self-energy is a single polarization bub¬ 
ble, fares in this system. In order to answer this, we 
start with the symmetric solution for which the contour 
plots of Fig. 5 show that both energy scales of the ex¬ 
act solution are reproduced for the interaction strengths 
considered. However, in the low energy scale, we only 
observe [P] and do not find a clear signature of a peak 
around uj/ujq ~ 1 for A ~ 2 for the propagation times 
accessed in this work. In the high energy scale, as the in¬ 
teraction is increased the fully self-consistent spectra be¬ 
come denser with non-uniformly separated peaks which 
do not move as a whole to higher energies. These obser¬ 
vations are all in a disagreement with the exact results 
which show uniformly two bare phonon frequency sepa¬ 
rated peaks moving to higher energies. This observation 
is however consistent with the previously discussed prop¬ 
erties of the approximate electron propagator for strong 
interactions. The asymmetric solution, once it is found, 
is observed to approach the non-interacting result i.e. the 
lowest frequency approaches the bare phonon frequency 
and higher lying structure looses intensity as the interac¬ 
tion is increased. This is expected since there is no room 
for particle-hole excitations in the symmetry-broken sys¬ 
tem, and thus the polarization bubble should tend to zero 
when the interaction is increased. As in the electronic 
case, also these spectra should fulfill sum rules given in 
terms of the zeroth and first moments by 

OO 

d<°> = - / tvD>(u) 

J 27 n 

— OO 

= trA, (13a) 

OO 

d (1 -* = — J wtr(o;D > (a;)) 

— OO 

= 2 E p q + E ep Q + tr(aJ7 J ^) , (13b) 

where A = tD M ( 0 + ), and E p q and E ep Q are the quan¬ 
tum contributions to the phonon and electron-phonon 
interaction energies defined in Ref. . The top panels 
of Fig. 5 show that these sum rules are approximately 
obeyed, and therefore an important consistency rela¬ 
tion is satisfied. Here it is noteworthy that although 
there is a clear change in the phonon energy, see the 
zeroth frequency moment, in the exact and fully self- 
consistent solutions, only the former displays a clear kink 
at A ~ 1.3 — 1.5 in the first frequency moment. Finally, 
the top panels of Fig. 5 highlight the lowest excitation 
energy labeled with [P] which is the dominant part of 
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FIG. 5. The exact (ED) and approximate (H, Gd, GD) 
phonon propagator as a function of the interaction A and fre¬ 
quency to. The top and bottom figures relate to the adiabatic 
ratio 7 = 1/2 and 7 = 1/4, respectively. The contour plots 
show |£>>(w)| /T on logarithmic scale for the symmetric (s) 
solutions and, starting from the critical interaction marked 
with a vertical dashed line and its value Ac, for the asym¬ 
metric (a) solutions. The middle panels show the zeroth d (() ' 
and first d (1 ' moments (points) of the spectra and the corre¬ 
sponding expectation values (lines) to illustrate fulfillment of 
Eqs. (13). The top panels show the intensities and positions 
of the lowest energy peak of |Dn(w)| /T for the symmetric 
(solid line) and asymmetric (dashed line) solutions labeled as 
[P] in the contour plots, (color online) 
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the spectra. The exact results show that this peak ap¬ 
proaches zero energy, but never actually reaches it, and 
gains intensity as a function of the interaction. This is 
in contrast to the non-interacting propagator in which 
this sole feature remains at the bare phonon frequency. 
The self-consistent Born approximation however captures 
both effects reasonably accurately up to A ~ 1.5 and gives 
a qualitatively similar trend even beyond it for the inter¬ 
actions considered here. 

To summarize, we found that the non-interacting prop¬ 
agator used in the Hartree and partially self-consistent 
Born approximations is an adequate approximation for 
the interacting phonon propagator only for very weak 
1 C 1 interactions. The fully self-consistent Born ap¬ 
proximation, on the other hand, captures the dominant 
low energy peak well up to borderline strong A ~ 1.5 
interactions and is therefore in a good agreement with 
the exact results in this range of interactions. It how¬ 
ever does not describe the second low energy excitation 
at wo, and therefore reproduces qualitatively only one of 
the signatures of a bipolaronic state observable in the 
phonon propagator. Finally, we remark that the absence 
of a peak at the bare phonon frequency for strong in¬ 
teractions is a likely factor for the observed too dense 
frequency structure of the electron propagator. 


B. Linear Response Functions 

The dynamics of a system of electrons and nuclei is in 
many cases dominantly determined by a linear response 
function provided that the system is perturbed suffi¬ 
ciently weakly. Many spectroscopic methods essentially 
rely on measuring these functions which makes them im¬ 
portant for understanding experiments. Here we investi¬ 
gate in particular the first order density-density response 
function of our model system. 


1. Method 


7 = 1/2 



-5 0 5 


Relative Displacement (u) 




Let us begin by explaining how we in practice calculate 
linear response functions in our time-dependent formal¬ 
ism. This is a prerequisite for understanding when it 
is reasonable to do linear response by time-propagation. 
The density-density response function is calculated by 
perturbing the system with the time-dependent poten¬ 
tial 

= ( 14 ) 

icr 

where Vi is the magnitude of the perturbation, and 6 is 
the Dirac delta function. Then we record the resulting 
spin-summed density 


FIG. 6 . The adiabatic potential energy surfaces Eo(u), E\(u ) 
and E 2 (u) for the three singlet eigenstates of the electronic 
clamped nuclei Hamiltonian as a function of the interaction 
A and relative displacement u, see Ref. 1 for details. The 
top and bottom figures correspond to the adiabatic ratios 
7 = 1/2 and 7 = 1/4, respectively. The left panels contain 
A = 0, 0.5,1.0,1.5, 2.0 cross-sections of the potential energy 
surfaces shown as contour plots in the right panels, (color 
online) 

which in the linear response regime satisfies 

Srii(t) = rij(i) - n-°V) 

3 


n i(t) = -2 


( 15 ) 
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where v is the norm of a vector with Vi as its compo¬ 
nents, nf\t) is the density of the unperturbed system, 
and the retarded component of 

the first order density-density response function. The 
response function is then given by 

, ( 16 ) 

v—0 

which we in practice evaluate by using the difference quo¬ 
tient ( rii(t ) — n[°\t))/vj with Vj sufficiently small and 
Vk = 0 for k j. Lastly, it is important to under¬ 
stand that applying the delta function potential amounts 
to choosing a new initial state after which the time- 
evolution is induced by the unperturbed Hamiltonian. In 
exact diagonalization, this is achieved by preparing the 
new initial state 

K)=e-^ v ^\*») , (17) 

where I’L^) is the N electron ground state, and which 
is subsequently propagated in the absence of the pertur¬ 
bation. In the Kadanoff-Baym equations, on the other 
hand, the same is achieved by choosing the new initial 
electron propagators 


G|(0;0) = e-^-^G^(0 ± ), 

(18) 

G} j (0;r) = e-^Gff(-r), 

(19) 


where Gfj(r) is the solution to the equilibrium Dyson 
equation. The electron and phonon propagators are then 
obtained by time-propagating the unperturbed Kadanoff- 
Baym equations. 

2. Stability 

The method described above is expected to work if 
the perturbation expansion of Eq. (15) is valid for the 
time scales of interest. It can however be that a possibly 
finite time-scale in which the expansion is good cannot be 
extended to cover the entire time-scale of interest. This 
can signal e.g. an unbounded linear response function. In 
the following, we show that this is the case for the Hartree 
approximation, and subsequently investigate whether or 
not the Born approximations show a similar behavior. 

The Hartree Eqs. (10) are a closed set of ordinary dif¬ 
ferential equations for the phonon field expectation value 
</>p(t) and the electron propagator In the two- 

site, two-electron Holstein model these equations can be 


rewritten as 

h = 4t kin T 2 , (20a) 

f i = -2 guT 2 , (20b) 

f 2 = -£kin« + 2pitTi, (20c) 

u = ujqP , (20d) 

p = —coqu + 2gn , (20e) 


where we have suppressed the time arguments and the 
overhead dot denotes the time-derivative. Moreover, 
n = n la - n 2a , u={u i - u 2 )/\/2, and p = (pi -p 2 )/V 2 
are the relative spin density, displacement, and momen¬ 
tum, while T i and T 2 are the real and imaginary parts of 
the density matrix element 712 = —iGf 2 , respectively. As 
shown explicitly in App. B, the density-density response 
function is the solution to the corresponding linearized 
equations of motion. If the linearization is performed 
with respect to an equilibrium solution which is a stable 
fixed-point, in the sense of Lyapunov ' 5 , of the original 
equations then the eigenvalues of the resulting Jacobian 
matrix have non-positive real parts . Moreover, if there 
are no repeated zero eigenvalues, then the zero solution 
of the linearized system is stable and furthermore any so¬ 
lution is bounded . In particular, the density response 
function is then bounded, that is 3 M > 0 independent 
of t such that |xf^(£)| < M for all t > 0. In order un¬ 
derstand when this is the case, we investigate below the 
stability of the fixed-points of the Hartree equations. The 
fixed-points whose stability is to be studied are just the 
symmetric 

n s = 0, 

r s , 2 = 0, 

u s = 0 , 

p s = 0 , ( 21 a) 

and asymmetric 

n a ± = ± \/l - A “ 2 , 

r a ±,2 = 0 , 

U a ± = 2 gn a ±/uj 0 , 

Pa± = 0. (21b) 

solutions of the equilibrium Hartree equations of 
Eqs. (A2) derived in Ref. 1. These equations are sub¬ 
ject to two constants of motion as both the eigenvalues 
of the reduced density matrix, which are either one or 
zero, and the total energy are conserved and give the 
constraints 

i = n 2 +A(r 2 1 + r 2 2 ), 

E= U j{p 2 +u 2 ^ 1 ) - dtkinTi - 2 gnu , ( 22 ) 

respectively. Then by following , and motivated by 
the first constraint, we introduce the coordinates 

z = 2Ti, 

n = \Jl — z 2 cos(d ), 

T 2 = \/l - z 2 sin( 0 )/ 2 , 

which represent cross-sections of the unit sphere with a 




drn(t) 

dvj 
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plane. The transformed equations of motion 

uz cos (9) 


9 — —2f kin + 2 g 


Vl — 


z = —2 gu\J 1 ~ z 2 sin(0), 

u = UJ 0 p , 

p = — ujqu + 2 < 7 \/1 — z 2 cos (9). 


and the total energy 

E = ^(p 2 + u 2 — l) — 2f kin ,z — 2gu\J\ — z 2 cos (9 ), 

then correspond to a Hamiltonian system with canoni¬ 
cal conjugate variables (0,z) and ( u,p ). The canonical 
transformation 


qi = —\/2(l - z) sin(0), 
pi = \/2(1 - z) cos (9), 

92 = ~P, 
p 2 = u , 

transforms this system into two non-linearly coupled os¬ 
cillators described by 

E= ^(pl + ql) + 4in (p\ + q\) 

- %gp2Pi \Jl - (pf + g?)/4, 

where we dropped a constant energy shift. This system 
of equations is a special case of the Hamiltonian system 
studied in which arises from the semi-classical equa¬ 
tions of the Dicke model ’ . Let us then denote xi = qi, 
X 2 = 92 , %3 = Pi, and X 4 = P 2 - The fixed points of this 
system are just related by coordinate transforms to the 
fixed points of Eqs. (21). At the symmetric fixed-point, 
we find the Hessian matrix 

/2f kin 0 0 0 \ 

S "o 2f k i n -2, 

\0 0 -2g oj 0 ) 

which is positive-definite for A < 1 and indefinite for 
A > 1, while at the asymmetric fixed-points, the Hessian 
matrix 


VVE a± 


/^kin(l + A) 0 
0 W 0 

0 0 

v 0 0 


0 

0 

4*kinA 

1+A- 1 
2\Z2gA — 1 
x/l+A-i 


0 

0 

2y/2gX~ 1 

x/l+A-i 

WO 


is positive-definite for A > 1. The symmetric equilibrium 
and asymmetric equilibria are then due to the Lagrange- 
Dirichlet theorem ’ 1 stable for A < 1 and A > 1, respec¬ 
tively. Moreover, since det(WA s ) < 0 for A > 1 also 
det(JVVi? s ) < 0, where J is the standard symplec- 
tic matrix and therefore the Jacobian matrix of the 


linearized Hamilton’s equations has an eigenvalue with 
a negative real-part. This implies that there also ex¬ 
ists an eigenvalue with a positive real part which means 
that the equations are linearly and nonlinearly unsta¬ 
ble . The symmetric equilibrium is therefore unstable 
for A > 1. The zero solution losing its stability while 
two new stable equilibria arise is a standard bifurcation 
known as the supercritical pitchfork bifurcation . The 
stability together with the fact that the Hessian matri¬ 
ces do not have zero eigenvalues for A ^ 1 implies that 
the response functions obtained for A < 1 and A > 1 us¬ 
ing the symmetric equilibrium and asymmetric equilibria 
are bounded functions. The linear instability of the sym¬ 
metric solution for A > 1 leads, on the other hand, to an 
unbounded response function as shown in App. B. 

This answers the question when it is in this context ap¬ 
propriate to do linear response properties at the mean- 
field level, but does not resolve this issue for the cor¬ 
related approximations. In this case, one cannot re¬ 
cast the equations as a set of ordinary differential equa¬ 
tions, but instead must consider the full two-time integro- 
differential equations which have non-linear integral ker¬ 
nels. As we are not aware of stability theory for such 
dynamical systems and it would go beyond the scope of 
the present work, we only resort to a working measure 
which is in the spirit of the stability of the equilibrium 
solutions. The working measure chosen here is a practical 
one: we introduce the norm 

IIHL = t max I n(t) - n(0)| , (23) 

compare it to the magnitude of the perturbation v, and 
if they remain in the same order of magnitude, we sug¬ 
gest that the equilibrium is stable. We emphasize that 
this measure is not equal to the stability even in the case 
of ordinary differential equations, but does give a practi¬ 
cal estimate whether or not a linear response calculation 
makes sense for the time scales accessed in this work. 

Having said this, Fig. 7 shows this norm for exact di- 
agonalization (ED) and many-body theory (H, Gd, GD). 
The results are obtained either by starting from a sym¬ 
metric, or asymmetric equilibrium solution. The exact 
results stay always in the same magnitude as the perturb¬ 
ing potential v = 10~ 3 which is expected since the time- 
dependent Schrodinger equation is linear. The other ex¬ 
treme is the symmetric mean-field solution for which the 
norm suddenly, although continuously as a function of 
the interaction, approaches one when the interaction ex¬ 
ceeds the corresponding critical value A = 1. At the same 
time, this norm remains in the same order of magnitude 
as the perturbation for the asymmetric mean-field solu¬ 
tion. This is consistent with the stability analysis pre¬ 
sented above illustrates that our working measure agrees 
in this case with Lyapunov stability. These cases give us 
some confidence in looking at the symmetric and asym¬ 
metric solutions of the Born approximations. The results 
for the asymmetric cases of these approximations indicate 
that they behave qualitatively similar to the asymmet¬ 
ric equilibria of the mean-field suggesting that they are 
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ED - H — Gd - GD — 




FIG. 7. The norm HdnH^ of Eq. (23) as a function 
of the interaction A for the perturbation of Eq. (14) with 
V\ = 10 _ 3 ,u 2 = 0. The top and bottom panels correspond 
to the adiabatic ratios 7 = 1/2 and 7 = 1/4, respectively. 
The solid and dashed lines correspond to the symmetric and 
asymmetric solutions, respectively. The vertical dashed lines 
and the associated values A c denote the critical interactions 
which are ordered according to to H, Gd, and GD from left 
to right, (color online) 


stable. This qualitative agreement remains true for the 
symmetric equilibria of the partially self-consistent case 
in which the norms approach one when passing the crit¬ 
ical value of interaction. The norm of the symmetric so¬ 
lution of the fully self-consistent approximation however 
remains in the same order of magnitude as the perturba¬ 
tion. We note that the declining norm for 7=1/4 and 
higher interactions is in the exact case related to the fact 
that the propagation length is shorter than the period of 
the lowest excitation of the system. 


These observations together with the interpretation 
given to our working measure of stability then suggest 
that the partially self-consistent approximation has the 
same qualitative stability properties as the mean-field. 
On the other hand, both equilibrium solutions of the fully 
self-consistent Born approximation are observed to be 
stable in the sense of our working measure. This means 
that, on the contrary to the Hartree and partially self- 
consistent Born approximation, it is possible to do lin¬ 
ear response with respect to the symmetric equilibrium 
of the fully self-consistent Born approximation for the 
time-scales addressed here. 


3. Bethe-Salpeter equation 


The Bethe-Salpeter equation is the standard approach 
to calculate linear response functions in many-body per¬ 
turbation theory . Here, we discuss the connection 
between the many-body approximations used in this 
frequency-domain approach and in the time-dependent 
approach applied in the present work. We start by noting 
that the density response function of Eq. (16) is the re¬ 
tarded = E av' Xfaiajv'ja'fa 0 ) component of the 

generalized, contour-time response function 


Xij,ki( z 'i z ') ~ ^ Tr 


Tie 


D -iJ c dzH(z 


) A‘j ij (z)A A / k i(z') 


(24) 


where 7 = c^Cj is the one-body reduced density matrix 
operator and A7^ the corresponding fluctuation oper¬ 
ator. Note that we switched here to collective indices 
containing both spatial and spin degrees of freedom. In 
the standard approach, the generalized response function 
satisfies the equation 


Xij,ki{ z i z') = Pij,ki(z; z') + ^ y / dzdz Pij, rs (z; z) 

PQ rstu 

x Mf r (z)dp Q (z]z)M^ u (z)xut,ki{^ z ') > 

(25a) 

Pijikii^z, z ) — t ^ ^ I dzdz Gip(z, z'j 

pq J c 

x Gqj{.z, z)T pq -ki{z,z, z ), (25b) 

where Pij k i{z\ z') is the irreducible polarizability defined 
in terms of the irreducible vertex function Tij^i{z, z'\ z"). 
These equations are valid for any many-body approxima¬ 
tion which includes the mean-field, Hartree term while 
beyond mean-field effects are incorporated into the irre¬ 
ducible vertex function. This function satisfies the Bethe- 
Salpeter equation 3 

r z'\z") = 5u5 jk 5{z , z')8(z', z") 

+ Jdzdz 1 dzdz 1 Ki q ^ p j(z, z\ z ', z') 

pqrs q 

X G pr {z'\z)G sq (z'\ Z^)^rs : kl (Z,Z ,Z ) , 


where the four-point integral kernel is defined as 


Kij,ki{z,z'\z,z') 


<5E X c,ii (-2; z ) 
5G k j{z; z') 


(26) 


with the subscript xc denoting the exchange-correlation 
self-energy. The diagrammatic form of this equation is 
shown in the top panel of Fig. 8. It has been shown that 
a density response function obtained by time-propagation 
of the Kadanoff-Baym equations with a self-energy £ is 
equivalent to a solution of Eqs.(25) with the vertex satis¬ 
fying the Bethe-Salpeter equation with a four-point ker¬ 
nel of Eq. (26). Thus by calculating the response function 
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FIG. 8. The Bethe-Salpeter equation for the irreducible ver¬ 
tex function (top) and functional forms of the four-point ker¬ 
nels (bottom) in the Hartree (H), and the partially (Gd) and 
fully (GD) self-consistent Born approximations. A line with 
an arrow indicates a dressed electron propagator, while single 
and two-fold wiggly lines represent bare and dressed phonon 
propagators, respectively. An open circle represents a connec¬ 
tion for a phonon propagator and a closed circle (kernel) or a 
dashed line (BSE) a connection for an electron propagator. 


via time-propagation using Eqs. (15) and (16) we arrive 
by current standards at a high-level solution of the Bethe- 
Salpeter equation. 

Figure. 8 shows diagrammatically the approximate 
four-point kernels related to the self-energy approxima¬ 
tions used in the present work. In the Hartree approxi¬ 
mation, the irreducible vertex function is the bare vertex, 
and the response function is hence the sum of all bub¬ 
ble diagrams. The partially self-consistent Born approx¬ 
imation leads to a vertex function which consists of all 
the ladder diagrams and to the bubbles-and-ladders se¬ 
ries for the response function. In addition to such terms, 
the fully self-consistent Born approximation contains also 
two second order kernel diagrams in terms of the phonon 
propagators. These higher-order terms are not routinely 
considered in the zero-frequency • nor fully frequency- 
dependent cases 10 . Lastly, we note that despite of the so¬ 
phistication of these approximations there are in general 
no guarantees of their superiority over the conventional 
approximations. It has also been shown, that similar ap¬ 
proximations in the purely electronic case, can lead to 
undesired features like non-positivity 82 . 


4. Density-Density Response Function 


We have presented a stability analysis in order to un¬ 
derstand when it makes sense to calculate linear response 
properties in the context of the present work. Addition¬ 
ally, we have discussed the diagrammatic meaning of the 
response function obtained in this manner. Here we focus 
on the numerical results, that is analyzing the density- 
density response function obtained by time-propagation 
of the Kadanoff-Baym equations. In order to carry-out 
this analysis, we start by considering the exact response 
function which has the frequency-domain Lehmann rep¬ 
resentation 


M = 

n 


h N h N - 

n m n n] 

w + **7 


hnjhnj \ 
w + Q.% + 117 ) 


where = E^ — Eq is a neutral excitation energy, and 

= (\E'o r | Yjv I^r^) the corresponding real-valued 
oscillator strength. This form is used to analyze the exact 
(ED) density response function which is shown together 
with the approximate (H, Gd, GD) response functions 
in Fig. 9. The main panels contain contour plots which 
illustrate the overall structure of the spectra. The exact 
results show that the non-interacting response function 


(_l)i-i/2 (— iy-i/2 

uj - 2 f kin + ir] uj + 2 f kin + ir] ’ 


which consists of a single peak for positive frequencies, 
develops as a function of the interaction to a function 
comprising multiple excitation energies. In the case of 
weak A < 1 interactions, these excitations can be rea¬ 
sonably well identified as phonon sidebands i.e. as multi¬ 
phonon excitations either from the non-interacting elec¬ 
tronic ground or first, singly-excited, excited state. The 
sidebands corresponding to either electronic state are 
separated roughly by two bare phonon frequencies for 
weak interactions. Moreover, there is an extremely faint 
peak located at w/£ k i n ~ 4.5 for 7 = 1/2 coinciding en¬ 
ergetically with the non-interacting doubly-excited elec¬ 
tronic state plus a single phonon. As the interaction 
is increased, the structure associated initially with the 
singly-excited electronic state labeled with [E\ in the fig¬ 
ures moves as a whole to higher energies. At the same 
time, the peaks related initially to the electronic ground 
state approach a bare phonon frequency separated dis¬ 
tribution with the lowest excitation labeled with [P] in 
the figures approaching zero energy and gaining inten¬ 
sity. These spectral features can be understood from the 
adiabatic potential energy surfaces of Fig. 6 as discussed 
in Sec. IV A 2. In the following, we instead focus on the 
dominant low energy peak [P] with the aim to identify 
it as a signature of a bipolaronic system. We start by 
writing the exact time-dependent density as 

/»oo 

rii(t) = du (2 Pu(u]t) + Pi 2 (u-,t)) , 
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7 = 1/2 

ED - H - Gd - GD - 



Interaction (A) 


7 = 1/4 

ED - H - Gd - GD - 



Interaction (A) 



Interaction (A) 



0 1 2 0 1 20 1 20 1 22 
Interaction (A) 


FIG. 9. The exact (ED) and approximate (H, Gd, GD) retarded density-density response function as a function of the 
interaction A and frequency w. The top and bottom panels correspond to the adiabatic ratios 7 = 1/2 and 7 = 1/4, respectively. 
The contour plots on the right show |trx fl (w)|/T for the stable symmetric (s) or asymmetric (a) solutions. The small top 
panels show the first n (1) moment (points) and the corresponding expectation value (lines) to illustrate fulfillment of Eq. (28). 
The left panels show the intensities and positions of the two peaks of |trx fl (w) | /T labeled with [P] and [E] in the contour 
plots. The stable solutions of H and Gd are denoted with a solid line irrespective of their symmetry. The stable symmetric and 
asymmetric solutions of GD are denoted with a solid and dashed line, respectively, (color online) 


where Pij(u;t) defined by 


Pu(u-,t) = \(i t,U;«|^ =2 (t)}| a , 
Pi2(u;t) = Y l( lo ',2cr , ;w|'i'^ =2 (t))| 2 , 

GO 


written as 



du Quj{u\t ), 


where 


efkjfat) 


dPik (u; t) 


dvj 


v—0 


(27) 


is the time-dependent joint probability to find the elec¬ 
trons at sites i and j and nuclei at the relative coordi¬ 
nate u at time t. Here we use the notation |\E/^ =2 (f)) = 
exp(-iH M t)\^g= 2 ) with |^^ =2 ) defined in Eq. (17) and 
\ia,ja';u ) = c\ a c^ a , |0) e |u) where |0) e is the electronic 
vacuum and [u) is an eigenstate of u. The exact den¬ 
sity response function can be then according to Eq. (16) 


is a response function describing how the ground state 
joint probability P i? (u: 0) changes as a function of time 
due to a weak perturbation. The response function 
012 does not contribute to the density response 

function since it is an odd function under the interchange 
u 77 —u as follows from the full inversion symmetry 
of the model. In Ref. 1, we use the ground state joint 
probabilities as ingredients of a working definition of a 
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dominantly bipolaronic ground state. The probabilities 
Pij(u ; 0) shown in the top panels of Fig. 10 illustrate the 
fact that as the interaction increases one is most likely to 
find the system in a state in which both electrons occupy 
the same site with an accompanying nuclear displace¬ 
ment. In particular, at A = 2.0 for 7 = 1/2 and A = 1.7 
for 7 = 1/4, the ground state of the system has according 
to Ref. 1 crossed over to a dominantly bipolaronic state. 
Next, we illustrate how these distributions behave in the 
linear response regime by showing the time-average 


( Qik,j)(u ) 



I Qik,j {u'i t) | J 


as a function of the interaction and displacement in the 
left contour plots of Fig. 10. The final time T is cho¬ 
sen here so that T/t ^ « 470 and T/t ^ « 9360 for 
7=1/2 and 7 = 1/4, respectively. The results indicate 
that i) 011,1 (u, t) and 022,1 (u, t) are on average larger than 
012,1 (u, £), and that for a sufficiently strong interactions 
the latter become suppressed while the former gain mag¬ 
nitude. The maxima max„,t e ro,T] |0ifc,j(tt;01 shown in 
the insets underneath the averages further support these 
statements. Moreover, we observe that when the ground 
state distributions become spatially polarized as the in¬ 
teraction is increased, also the response functions follow 
the same trend. Thus we find that ii) the spatial shapes 
of the initial distributions remain qualitatively invariant 
in the linear response regime as a function of time. In 
order to illustrate the temporal behavior of the dominant 
response functions 0u,i(w;t), we show them in the right 
panels of Fig. 10 for the initially bipolaronic systems at 
A = 2.0 and A = 1.7 for 7 = 1/2 and 7 = 1/4, respec¬ 
tively. Firstly, these results agree with the conclusion 

ii) on the spatial structure, and secondly, they show that 

iii) probability density is redistributed between 0n,i(w; t) 
and 022,1 (u; t) mainly on a time-scale given by the energy 
scale of [P], while the energy scale given by [E] is seen as 
superimposed small amplitude oscillations. The points i), 
ii) and ii) combined allow us to conclude that, in agree¬ 
ment with the working definition of Ref. !, the system is 
in a dominantly bipolaronic state at each instant of time. 
Moreover, we understand the oscillation of the probabil¬ 
ity density between 0n,i(u;t) and 022,1 (u;t) to represent 
the motion of a bipolaron appearing according to iii) on 
a time scale set by [P]. Finally, this is seen in the density 
response function according to Eq. (27) as the emergence 
of the dominant low energy excitation [P]. 

As we have described the exact density response func¬ 
tion, we are prepared to investigate how the many-body 
approximations describe it in order to understand their 
limitations. Let us begin with the mean-field, Hartree ap¬ 
proximation in which the density response function can 
be obtained analytically as the solution to the linearized 
Hartree equations as shown in App. B. This response 
function is given by 
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FIG. 10. The exact (ED) joint-probabilities Pq(w; 0) and re¬ 
sponse functions Qik,j{ w; t) shown for the adiabatic ratios 7 = 
1/2 and 7 = 1/4 in the top and bottom figures, respectively. 
In the top panels, we show P,;, (u; 0) as a function of the dis¬ 
placement u for the interactions A = 0.2,1.7, 2.0 whose color 
code is shown in the bottom panels. The middle panels con¬ 
tain (Qij,i) as a function of the displacement and interaction, 
and the bottom panels display max u, t £ [0, T] | 0^,1 (u; t )| as a 
function of the interaction. In right panels, we show 0 ii,i (u; t) 
as a function of the displacement and time t for the interac¬ 
tions A = 2.0,1.7 with regions closed by red (blue) contour 
lines being positive (negative) values, (color online) 
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where 


X± = 


2tkin(wo - W± 2 )' 


/ .IL 

UJ ± 


(w 2 -u 4 2 ) +4A w 2 t 2 in 


are the oscillator strengths and 


<4 = 


N 


4/ 2 

^ 6 kin 



/ ie^t 2 in (A-Ty \ 

v (-o 2 +4tD 2 y 


w ± = 



/ 16 a,gtg in (A^ ^ry \ 

V K+4iLA 2 ) 2 J’ 


are the frequencies for the symmetric (s) A < 1 and asym¬ 
metric (a) A > 1 solutions, respectively. The Hartree 
response functions shown in the subpanels H of Fig. 9 
thus consists of two contributions: the high [E] (x+ , w +) 
and low [P] (x",w") energy peaks related, for weak in¬ 
teractions, to the non-interacting electronic excited and 
ground states plus zero and one phonon, respectively. 
Firstly, we observe that as the interaction is increased, 
the initial distribution given by [E\ remains nearly invari¬ 
ant up to the critical interaction A = 1 beyond which its 
energy increases linearly as a function of the interaction. 
Secondly, the low energy peak [P], which has no inten¬ 
sity in the non-interacting A = 0 case, gains intensity 
and approaches the zero energy as the interaction is in¬ 
creased. Moreover, when it reaches the zero energy at the 
critical interaction A = 1, its intensity given by the oscil¬ 
lator strength diverges. As explained in Sec.IVB2, 
by increasing the interaction beyond this point, we make 
the symmetric equilibrium solution unstable, and there¬ 
fore we change the initial state to one of the asymmetric 
solutions. This leads again to a well-defined first order 
response with the intensity of the low energy peak becom¬ 
ing finite and its frequency approaching the bare phonon 
frequency as the interaction is increased. We understand 
these results in terms of the adiabatic potential energy 
surfaces so that the mean-field approximation captures 
the lowest adiabatic potential energy surface of Fig. 6 
becoming more shallow which leads to the lowest excita¬ 
tion approaching the zero energy. This agrees with the 
observation that, as the Hessian matrices of Sec. IVB 2 
indicate, there is a direction in the energy landscape in 
the neighborhood of the symmetric equilibrium solution 
such that along it, as A —> 1, the approximately har¬ 
monic energy surface becomes more shallow. Moreover, 
at A = 1, this harmonic surface becomes completely flat, 
and as the linearized equations describe only this local 
neighborhood, it appears as if exciting the system costs 
no energy which manifests itself as the divergence of the 
zero-frequency component of the response function. At 
this point, the lowest adiabatic potential energy surface 
forms the double-well structure which has for A = 1 also 
a locally flat energy landscape at u = 0 where its sec¬ 
ond derivative vanishes. Lastly, one can show that the 
Hartree ground state energy function is equivalent to the 


lowest adiabatic potential energy surface Eq(u ) by en¬ 
forcing in Eq. (22) that 2 gn = u and that Ti satisfies the 
equilibrium Hartree equations derived in Ref. 1. This 
suggests that the Hartree approximation captures the 
formation of the double-well potential but needs to fall 
into one of the two minima in order to minimize the en¬ 
ergy. By doing so, it sees again a nearly harmonic surface, 
which appears in Fig. 9 as the lowest excitation becoming 
finite and approaching the bare phonon frequency. 


Let us then discuss the density response functions ob¬ 
tained for the stable equilibrium solutions of the par¬ 
tially (Gd) and fully (GD) self-consistent Born approx¬ 
imations shown in Fig. 9. The partially self-consistent 
results shown in the subpanel Gd of Fig. 9 indicate that 
the main qualitative difference to the Hartree approxi¬ 
mation is that there is a sideband structure related to 
the excitation [E\. The sidebands are separated roughly 
by two bare phonon frequencies in agreement with the 
exact results for weak interactions, but do not move 
to higher energies as a function of the interaction as 
clearly as the exact spectra does. Moreover, when A 
exceeds Ac, the ground state becomes asymmetric and 
new symmetry-forbidden excitations emerge in-between 
the original sidebands. In the low energy scale, we in¬ 
stead do not observe new qualitative differences to the 
mean-field solution, in particular we still find that at Ac, 
the low energy peak [P] reaches the zero energy with its 
intesity diverging. The symmetric solution of the fully 
self-consistent Born approximation does, however, show 
a qualitative difference as shown in the subpanel GD (s) 
of Fig. 9. In the low energy scale, we observe that as the 
interaction is increased, the low energy peak [P] moves 
initially towards the zero energy but, in contrast to the 
mean-field and partially self-consistent results, does not 
reach it for the parameters considered in this work. This 
is in an agreement with the exact solution in which, how¬ 
ever, the lowest excitation becomes increasingly close to 
the ground state, while in the fully self-consistent approx¬ 
imation, we observe that it approaches a finite non-zero 
value. In the high energy scale, we on the other hand ob¬ 
serve that as the interaction is increased, the peaks above 
[E\ become non-uniformly spaced and too dense in com¬ 
parison to the exact solution. These shortcomings, as 
well as the fact that the spectra do not move apprecia¬ 
bly to higher energies as a function of the interaction, 
are similar to what we observed for the equilibrium elec¬ 
tron propagators in Sec. IV A 1. Finally, the asymmetric 
solutions shown in the subpanels GD (a) are similar to 
the partially self-consistent solutions for A > Ac except 
for the additional excitation at u>/u >o ~ 2. The low en¬ 
ergy sidebands seen in the exact solution are then likely 
merely symmetry-forbidden in the symmetric solution of 
the fully self-consistent approximation. As the last re¬ 
mark on the overall structure, as shown in the top insets 
of Fig. 9, the development of these spectra as a function 
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of the interaction is consistent with the f-sum rule 


n« ee - 


OO 

/ duj o, . 

— wtrx (w) 

m 


= -2 (Ee - El° c ) - 2 (E ep - E%) , (28) 


which relates the first moment of the density-density re¬ 
sponse function to the electron ( E e ,E l ° c ) and electron- 
phonon interaction energies (E ep ,E l ° p ), where super¬ 
script loc refers to the site-diagonal part of the corre¬ 
sponding energy. The fact that also the many-body ap¬ 
proximations satisfy this sum rule is proven for the purely 
electronic case in • . 

Next, we focus on the two most significant features of 
the density response function for the weak and interme¬ 
diate interactions. These are the intensity and position 
of [P] and [F] shown in the left panels of Fig. 9. In the 
exact case, as the interaction is increased, \E\ loses mag¬ 
nitude and moves towards higher energies while [P] gains 
intensity and approaches the zero energy. The former is 
dominant up to borderline strong A ~ 1.5 interactions 
although the latter is appreciable already for A ~ 1, and 
its impact to the response properties is emphasized by its 
different energy scale. In the Hartree approximation, we 
observe that, as a function of the interaction, [E\ is nearly 
invariant for A < 1 and that [P] behaves in a divergent 
manner as desribed above. The mean-field approxima¬ 
tion therefore agrees with the exact solution only for very 
weak interactions A C 1. The partially self-consistent 
Born results are qualitatively similar to the mean-field 
results. On a quantitative level, it reproduces the exact 
intensity and position of [E\ better but deviates consid¬ 
erably for intermediate A ~ 1 interactions. This together 
with the observed divergence of [P] implies that it can 
be said to agree well with exact resutls only for weak 
A < 1 interactions. In both approximations, we find that 
the intensities of [F] and [P] decrease rapidly as a func¬ 
tion of the interaction once A exceeds Ac- The density 
response to a weak parturbation is thus suppressed for 
a sufficiently strong interaction which is consistent with 
the localized nature of the asymmetric equilibrium solu¬ 
tions discussed in Ref. . Moreover, we find that also 
the asymmetric solution of the fully self-consistent Born 
approximation behaves qualitatively in this manner for 
the interactions it has been found. Lastly, the symmetric 
solution of the fully self-consistent Born approximation 
reproduces the exact positions and intensities of [P] and 
[P] well up to intermediate interactions A ~ 1 with the 
intensity of [P] being good also for stronger interactions. 
In this approximation, the low energy excitation [P] does 
not reach the zero energy nor does it diverge as a func¬ 
tion of the interaction, but its position and intensity do 
not still agree even qualitatively with the exact solution 
for strong A > 1 interactions. 

To summarize, we have found similarly as in 
Sec. IV A 1 that the Hartree, and partially and fully self- 
consistent Born approximations are in a good agreement 
with the exact results for very weak A -C 1, and weak 


A < 1 and up to intermediate A ~ 1 interactions. In par¬ 
ticular, we have shown that the exact density response 
function has, for sufficiently strong A > 1 interactions, 
a dominant low energy excitation which none of the ap¬ 
proximation describe qualitatively correctly. Moreover, 
we have related this excitation to the response of a bipo- 
laron to a weak perturbation by analyzing it using the 
time-dependent joint probabilities. Instead of describing 
the low energy excitation, the Hartree and partially self- 
consistent Born approximations give rise to a divergence 
of the response function at the critical point Ac- This 
has been explained by relating it to the formation of the 
double-well structure in the lowest adiabatic potential en¬ 
ergy surface which we have also associated with the bipo- 
laronic crossover in Ref. . Finally, we have shown that 
the divergence can be prevented by dressing the phonon 
propagator self-consistently at the level of the fully self- 
consistent Born approximation. 


5. Phonon Propagator Revisited 

The density response function which we calculated 
and discussed above describes electron density fluctua¬ 
tions which in turn couple to nuclear density fluctuations 
described by the phonon propagator. This is formally 
shown by the Dyson equation of Eq. (4b) which can be 
written as 

D(z ; z') = d(z ; z') 

+ Jdzdz' d(z ; z)TL r {z\ z')d.{z'\ z '), (29) 

G 

n t,pq{z; z') = M^(z)xij,ki( z 'i z')M^(z '), (30) 

ijkl 


where the reducible self-energy n r is determined by 
the generalized response function of Eq. (24). This 
means that a density response function obtained by time- 
propagation can be also used to obtain a new phonon 
propagator. Here we use this relation to identify some 
phonon self-energies and discuss whether or not they lead 
to better nuclear properties than the fully self-consistent 
Born approximation (GD). In order to do this, due to 
computational reasons instead of using the equilibrium 
frequency-domain version of Eq. 29, we perturb the sys¬ 
tem with the instantaneous force 


F P (t) = 5(t)F P , 

where Fp is the magnitude of the perturbation. We then 
record the resulting phonon field expectation value which 
in the linear response regime satisfies 

54> P {t) = 4> P (t) - <j>p\t ) 

= 5>£ Q (t)F Q + 0(F 2 ), 

Q 
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where <j)p^ ( t ) is the expectation value of the unperturbed 
system, and Dpq(t ) is the retarded phonon propaga¬ 
tor. In the Kadanoff-Baym equations this perturbation 
amounts to choosing 

c/) P (0) = 4>p - i'^cypqFq , 

Q 


where is the equilibrium expectation value, as the 
new initial condition and subsequently solving the equa¬ 
tions of motion in the absence of this perturbation. The 
phonon propagator is then given by 


D$ Q (t) 


d<j>p{t) 


dF Q 


F =0 


which we in practice evaluate by using the difference quo¬ 
tient (i j>p(t) — (f>p\t))/F q with sufficiently small Fq and 
Fp = 0 for R ^ Q. It can be shown that this propagator 
satisfies Eq. (29) and its irreducible version of Eq. (4b) 
with the irreducible self-energy 


n pq(z-, z') = Mp(z)Pi jt ki{z\ z')M^(z '), 

ijkl 


where P is the irreducible polarizability of Eq. (25b). 
The phonon propagators obtained by time-propagation 
are thus related to irreducible self-energy functionals 
whose lowest-order diagrammatic expansions are shown 
in Fig. 11. The phonon propagator obtained in this man¬ 
ner in the Hartree (H), and partially (Gd) and fully 
(GD) self-consistent Born approximations are respec¬ 
tively given by 

n t d-H(z;z') = nB[GH,d](2;2 / ), 

n t d-Gd(^;2 / ) =n B h[GGd,d\(z-,z'), 

n t d-GD(z; z') = IIblx[Ggd, D GI )](z; z '), 

where we have introduced the prefix ’td-’ referring ot 
time-dependent to distinguish these self-energies from 
their original counterparts. This shows explicitly that 
these approximations are not self-consistent i.e. the prop¬ 
agator satisfying the Dyson equation is not the same as 
the ones in the self-energy diagrams. 

The results can be anticipated by noting that the non¬ 
interacting phonon propagator is a function peaked at 
the bare phonon frequencies. The reducible frequency- 
domain Dyson equation then suggests that the phonon 
propagator satisfying it has a similar frequency content as 
the density response function with weight redistributed 
around the bare phonon frequencies. This already gives 
a picture how good are the phonon propagators obtained 
from the density response functions of Sec. IV B 4. How¬ 
ever, let us try to make this picture more quantitative. 
Figure 12 shows the Fourier transforms of the retarded 
phonon propagators obtained by time-propagation for 
the many-body approximations (td-H, td-Gd, td-GD). 
The exact (ED) and fully self-consistent Born (GD) equi¬ 
librium propagators are also shown for reference. We 
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FIG. 11. The phonon self-energies, or their functional forms, 
corresponding to the phonon propagators obtained by time- 
propagation. The Hartree (H), and partially (Gd) and fully 
(GD) self-consistent Born approximations relate to the bubble 
(. B ), and bubble-ladders (BL) and bubble-ladders-exchange 
(BLX) self-energy functionals, respectively. A line with an 
arrow indicates a dressed electron propagator, while single 
and two-fold wiggly lines represent bare and dressed phonon 
propagators, respectively. An open circle represents a connec¬ 
tion for a phonon propagator. 


have discussed the reference results and their physical 
content in Sec. IV A 2 so here we focus directly to com¬ 
paring the different approximations. The contour plots 
show that none of the new approximations improve the 
qualitative description of the low energy peak [P] which 
dominates the spectra. Moreover, only td-Gd with a 
symmetry-broken ground state for A > Ac produces a 
sideband structure for this excitation. We also observe 
that td-GD has a slightly larger sideband separation for 
the electronic excitation when compared to the fully self- 
consistent Born (GD) spectra. Here we remark that the 
numerical results for td-H agree with the analytical re¬ 
sults and discussion on the phonon vacuum instability 
presented in Ref. 1. Lastly, we show the position and 
intensity of [P] relative to its exact position and inten¬ 
sity in the top panels of Fig. 12 for weak interactions 
A < 0.5. The results highlight, as expected in a pertur¬ 
bative regime, that td-H deviates the most and td-GD 
the least from the exact result. Furthermore td-Gd and 
GD give similar results in this regime with the former 
being slightly better as it includes all the diagrams up to 
fourth order in the electron-phonon interaction. 
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Overall due to the poor description of the lowest exci¬ 
tation, the td-H, td-Gd, and td-GD approximations are 
only valid in the regime of weak interactions A < 1 in 
which td-Gd and td-GD improve on GD. The qualita¬ 
tive behavior for larger interactions A ~ 1 shows that 
although one can obtain sophisticated many-body self¬ 
energies by means of time-propagation, they do not nec¬ 
essarily improve the description of the physics. In par¬ 
ticular, we find that infinite summation schemes for self¬ 
energy diagrams, which include vertex corrections, lead 
to deterioration of the spectral properties of the equilib¬ 
rium propagator evaluated with a single dressed polar¬ 
ization bubble for intermediate to high interactions. 


V. CONCLUSIONS AND OUTLOOK 

We have introduced a method based on time- 
dependent many-body perturbation theory aimed at 
studying interacting electrons and phonons. The many- 
body approximations used here are the Hartree (H), and 
the partially (Gd) and fully (GD) self-consistent Born 
approximations. The method has been applied to in¬ 
vestigate both the non-neutral and neutral excitation 
spectra of a two-site, two-electron Holstein model. We 
have presented results for the frequency-domain ground- 
state electron and phonon propagators, as well as for the 
density-density and displacement-displacement linear re¬ 
sponse functions. The results have been compared with 
numerically exact results obtained by exact diagonaliza- 
tion (ED) in order to assess their quality and relate a 
physical picture to the behavior of the many-body ap¬ 
proximations. 

In Ref. 1, we found that the approximations stud¬ 
ied here support multiple, concurrently co-existing solu¬ 
tions some of which exhibit a broken reflection symmetry. 
The asymmetric solutions were found once the electron- 
phonon interaction A reached a critical value Ac, and 
were understood to mimic the bipolaronic crossover of 
the exact solution. The asymmetric solutions were fur¬ 
thermore found to be the lowest energy solutions for a 
large range of parameters. The total energies, and nat¬ 
ural occupation numbers, also suggested that the sym¬ 
metric solution of the fully self-consistent Born approxi¬ 
mation describes partially the crossover to a bipolaronic 
state. In the present work, we studied the frequency 
structure of the ground state propagators for a restricted 
range of adiabatic ratios 7 = 1/2,1/4 and interactions 
As [0,2] which allows us to complete some of the obser¬ 
vations made in Ref. 1. Firstly, the frequency-integrated 
observables obtained from the electron propagator are in 
a better qualitative agreement with exact results than 
the frequency-resolved objects themselves. In particular, 
our results show that none of the approximations give an 
electron propagator in which there is a rigid shift or redis¬ 
tribution of the spectral weight comparable to the exact 
solution for high A > 1 interactions. The phonon propa¬ 
gator does not moreover show a clear spectral fingerprint 
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FIG. 12. The phonon propagators as a function of the in¬ 
teraction A and frequency w. The top and bottom figures 
correspond to the adiabatic ratios 7 = 1/2 and 7 = 1/4, re¬ 
spectively. The contour plots show Dii(u;)| /T for the exact 
(ED), denoted with red dots, and approximate (GD, td-H, 
td-Gd, td-GD) solutions. The top panels show the difference 
between the approximate and exact intensity and position of 
the lowest energy peak of |Df 1 (u;)| /T labeled with [P] in the 
contour plots, (color online) 
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of a double-well structure in the fully self-consistent Born 
approximation for the parameters considered. As all of 
these features have been identified as spectral signatures 
of a bipolaronic state, our results here favor the state¬ 
ment that none of the approximations describe even par¬ 
tially the bipolaronic crossover for the parameters con¬ 
sidered. This said, the results for the electron prop¬ 
agator can be roughly summarized by concluding that 
the Hartree, partially self-consistent Born, and fully self- 
consistent Born approximations agree with the exact re¬ 
sults up to very weak A <C 1, weak A < 1, and weak 
to intermediate A ~ 1 interactions, respectively. The 
non-interacting phonon propagator used in the Hartree 
and partially self-consistent Born approximations is only 
valid for very weak A< 1 interactions while the dressed 
propagator of the fully self-consistent Born approxima¬ 
tion agrees reasonably well with the exact results up to 
borderline strong A ~ 1.5 interactions. 

The linear response functions studied in the present 
work are obtained by time-propagation starting either 
from a symmetric or an asymmetric ground state solu¬ 
tion. In the case A > A c, these solutions co-exits and 
we find that the symmetric solutions of the Hartree and 
partially self-consistent Born approximations are unsta¬ 
ble, while the asymmetric solutions are stable, against 
a small asymmetric perturbation. The symmetric and 
asymmetric solutions of the fully self-consistent Born ap¬ 
proximation are, on the other hand, shown to be stable 
against the same perturbation. By identifying the sta¬ 
ble ground state solutions, we have been able to evalu¬ 
ate linear response functions corresponding to these so¬ 
lutions. In particular, the density-density response func¬ 
tion obtained in the Hartree and partially self-consistent 
Born approximations is shown to have a zero-frequency 
component which appears and its intensity diverges as 
A approaches Ac- In the Hartree approximation, this is 
caused by the build-up of a double-well structure in its 
ground state energy surface exactly at A = 1 and by the 
need of the mean-field to minimize the total energy. Our 
results further show that the fully self-consistent Born 
approximation does not have a similar divergence. The 
comparison of the exact and approximate density-density 
response functions confirms that the range of validity of 
the many-body approximations is roughly the same as 
for the case of equilibrium propagators. In particular, 
none of the many-body approximations were able to de¬ 
scribe the lowest excitation of the exact response function 
for strong interactions A > 1 for which it is the domi¬ 
nant feature of the exact response function. By analyz¬ 
ing this excitation, we further identified it as a signature 
of a bipolaronic system, and hence in agreement with 
the conclusions based on the equilibrium propagators, 
this suggests that the approximations do not describe the 
crossover to the bipolaronic state. Finally, we could by 
time-propagation obtain another phonon propagator as¬ 
sociated with a highly sophisticated self-energy, although 
a non-selfconsistent one. The results show that although 
the propagators obtained in the partially and fully self- 


consistent approximations are better than the equilib¬ 
rium propagator of the fully self-consistent Born approxi¬ 
mation for the perturbative Ad interactions, they lead 
to deterioration of the qualitative spectral properties for 
higher interactions. This suggests that, at least in this 
case, it is either important to maintain self-consistency, 
or that instead of infinite partial summations of dressed 
self-energy diagrams it is better to consider truncated 
approximations. 

In order to go beyond the approximations studied here 
in many-body perturbation theory one needs to intro¬ 
duce vertex corrections which are in particular needed to 
improve the properties of the electron propagator. This 
is realizable in an equilibrium theory but leads to a sub¬ 
stantial increase in the complexity of the time-domain 
method which makes the inclusion of vertex corrections 
challenging with the current computational resources. 
The two-time equations can be however cast as one-time 
equations via the GKBA ’ ’ which is a possi¬ 
ble way to overcome the numerical challenge and reduce 
more advanced approximations tractable. On the other 
hand e.g. in cavity quantum electrodynamics already our 
weak interactions are considered strong and thus the ap¬ 
proximations used here could be a valuable asset for in¬ 
vestigating non-linear time-dependent phenomena. In 
this context our results provide a basis for studying the 
approximation in an explicitly time-dependent situation, 
and for understanding properties of further approxima¬ 
tions e.g. the GKBA which could be used in order to 
address physics of larger or more realistic systems. 
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Appendix A: Hartree, TDSCF and Ehrenfest 

The Hartree approximation is shown here to be equiv¬ 
alent to the time-dependent self-consistent field (TD¬ 
SCF) approach and to the semi-classical Ehrenfest ap¬ 
proximation. This equivalence is shown by deriving 
the time-dependent self-consistent field equations, noting 
that they reduce to the Ehrenfest equations of motion, 
and finally by showing that their solutions can be used 
to construct the phonon held expectation value, and the 
electron propagator of the Hartree approximation. We 
begin by introducing the product ansatz 

\m) = \m) ix(*)>, 

where \4’(t)) and |y(f)) consist of only electronic and 
phononic degrees of freedom, respectively. By substitut- 
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ing this ansatz to the time-dependent Schrodinger equa¬ 
tion 

zd t |*(i)> = H(t) |*(i)> 

and projecting it to the states |x(f)) and |we arrive 
at the time-dependent self-consistent field equations 

id t | ip{t)) = Y ( hij{t) + Y M ij( t ) < t > p ( t )) eta, | tp(t)) , 

idt | x{t)) = [ ^ftpQ(i)<M i )<M*) 

V PQ 

+ YFp(t)4> P {t) + YK P j (*h ji(t)4>p) IxM) 

P ijP / 

where 

<pp(t) = (x{t )| ftp IxW). 

7 ij{t) = c)ci | ip(t)) , 


such that | ip(t)) can be written as a Slater determinant 
of the time-dependent orbitals ipi(t). In order to relate 
these orbitals to the electron propagator in the Hartree 
approximation, we further write down the equilibrium 


Hartree equations 


, 

(A2a) 

h% = h M + Y M P (t>P , 

(A2b) 

P 



(A2c) 


ij 


which are introduced in Ref. 1, and correspond to a set 
of non-linear eigenvalue equations for the eigenvalues tjf 
and eigenvectors ipjf . The electron propagator can be 
then written in the Hartree approximation in terms of 
the time-dependent orbitals ipi(t) obtained by solving 
Eqs. (Al) with cp 1 ^ and ipf 1 as their initial conditions. 
That is, in the Hartree approximation, the phonon field 
expectation values satisfy Eq. (Alb), and the electron 
propagator is given by 


are the phonon field expectation value and the reduced 
density matrix. In the derivation, we adopted the phase 
conventions 

IV>Cf)} = e^o dt ' (£p(t')-<x(t')b9p*(0>) # 

I x (t)) = e l ti 0 dt ' 5 

where we defined 

E e(t) = Y bijit) <^(i)l c\cj 1 1p(t)) , 
ij 

E p {t ) ^(*)l IxW) ■ 

PQ 

as the electron and phonon energies, respectively. The 
semi-classical Ehrenfest equations are then derived 
e.g. by introducing a polar expansion of the nuclear state 
and taking its classical limit . In our case however, the 
Heisenberg equations of motion for the phonon field ex¬ 
pectation values are equal to the classical equations of 
motion. Taking advantage of this property and the bilin¬ 
earity of the equation for | ip(t)} in terms of the electronic 
operators, we arrive at the Ehrenfest equations of motion 


idti>ij(t) = Y ( bjkit) 

k \ 

+ Y M j*i t ) < t >P i t ) S J' l l’ ik i t )’ ( Al& ) 

iY, a PQdtpQ(t) = Y J Cl P Q(t)4>Q(t) 

Q Q 

+ F P (t) + Y M Z(thj*(t)’ (^b) 

ij 


G£(i ;<) = - Y f+ifck^ljit'Wkiit) , 

k 

Gfjit-J) = -~Y f+ifck^tjit'Wkiit) , 

* k 

where /+ = 1 — /+, as readily verified by using Eqs. (Al) 
to check that Eqs. (5) with the Hartree self-energy are 
satisfied, and also by verifying that the Kubo-Martin- 
Schwinger boundary conditions J are met. 

Appendix B: Hartree Density Response Function 

Here, we calculate the density response function of 
the two-site, two-electron Holstein model in the Hartree 
approximation by applying the method discussed in 
Sec. IV B 1. In order to do this, we first rewrite the 
Hartree equations in a more convenient form by using 
the conserved total energy to reduce the number of de¬ 
pendent variables. That is, the total energy of Eq. (22) 
allows us to eliminate T 1; and subsequently by defining 
the vector 



we can rewrite the Hartree equations given in Eqs. (20) 
as 

x = f(x) 

( 4tkinX 2 
-ikin^i + 2ffar 3 ri(ari,a:3,a;4) 

U 0 X4 

-ujqX 3 + 2gxi 
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where 8t kin ri(:Ei, x 3 , x 4 ) = w 0 {x\ + x\ - l) - 45x10:3 - 
2 £'o- Here Eq denotes the total energy at t = 0 which 
is determined by the initial condition x°. The density- 
density response function of Eq. ( 16 ), which is in this 
case given by 


*£(*) = (- 1 ) 


i+i dn (t) 
dvj 


(B 2 ) 


i )—0 


is then the 1 = 1 component of the more general response 
function dxi/dvj |„=o which satisfies 


d dx 
dt dvj 


= J 


v—0 


dx 

dv.j 


v—0 


as seen by differentiating Eq. (Bl) with respect to vj. 
The Jacobian matrix Jij = d x . fi(x)\ v =o is a function of 
the unperturbed solution x\ v= o which is in our case ob¬ 
tained by propagating either the symmetric or asymmet¬ 
ric ground state solution of the equilibrium Hartree equa¬ 
tions given in Eq. ( 21 a) and Eq. ( 21 b), respectively. As 
these solutions are also fixed-points of Eq. (Bl), they are 
constant in time, and thus the Jacobian matrices for the 
symmetric (s) and both asymmetric (a) solutions given 
by 


Js 


J 


a 


0 

41 kin 

0 

0 

tkin 

0 

9 

0 

0 

0 

0 

w 0 

2 9 

0 - 

-w 0 

0 

0 

4 f kill 

0 


tkin 

A 2 0 

ff A ” 

-1 

0 

0 

0 


2 9 

0 

—Wo 



are time-independent. The initial conditions 

dx^/dvj | w= o for the symmetric (rj = s ) and asym¬ 
metric (77 = a) cases can be deduced from Eqs. ( 18 ) to 
be 


cobian matrix. The decomposition exists since the eigen¬ 
values of the Jacobian matrix given by zcu±, where 


\ 


w; 


+ 4 tl 


kin 


1 ± 


16w o i ki„( A - !) 


Wn 


' Kin) 


^n 2 + 4 tLA 2 


\ 


1 ± , 1 - 


16 ^ in ( A 2 -l) 

K+4iki„ A2 ) 2 , 


are non-degenerate for A 7! 1 . The diagonalizing similar¬ 
ity transformation given by 


X" = 


w^/ 4 t kil 

T) 

x + 


0- /4tki] 


-*w^/4t k i, 


rj ri 
—10J+X+ 




where x± = 2g/(l — w± 2 ), then allows us to write the 
response function as 


dx v {t) 

dvj 


v—0 


= Xt 


*X- 


05% 

dvj 


v—0 


where u> T) = diag (w+, wl , —w+, —wl) denotes a diagonal 
matrix. In particular, its first component according to 
Eq. (B 2 ) gives the density response function 


*£,,«(*) = -(-1 sin(a ,\t) + xl sin 


r) _ 

x± = 


2tkin(wo ^^± 2 Y 


V 

<4 


( uj 2 - u 1 2 ) +4A w 2 t 2 in 


where we introduced the Heaviside function to enforce 
the correct causal structure. 


dx° s 


dV;j 


V—0 


dx°a 

diu 


v—0 



-(Sij - ?2i)A-V4 


respectively. The equation for the response function is 
linear, and hence admits the solution 


dx v (t) 

dvj 


v—0 


= e 


dx® 
dvj 


v—0 


where we restored the explicit time-dependence. The 
task is then to evaluate the matrix exponential which 
is done here by using the eigendecomposition of the Ja- 
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